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The method of choice for describing attractive quantum systems is Hartree-Fock-Bogoliubov (HFB) 
theory. This is a nonlinear model which allows for the description of pairing effects, the main 
explanation for the superconductivity of certain materials at very low temperature. 

This paper is the first study of Hartree-Fock-Bogoliubov theory from the point of view of nu- 
merical analysis. We start by discussing its proper discretization and then analyze the convergence 
of the simple fixed point (Roothaan) algorithm. Following works by Cances, Le Bris and Levitt for 
electrons in atoms and molecules, we show that this algorithm either converges to a solution of the 
equation, or oscillates between two states, none of them being a solution to the HFB equations. 
We also adapt the Optimal Damping Algorithm of Cances and Le Bris to the HFB setting and we 
analyze it. 

The last part of the paper is devoted to numerical experiments. We consider a purely gravita- 
tional system and numerically discover that pairing always occurs. We then examine a simplified 
model for nuclcons, with an effective interaction similar to what is often used in nuclear physics. 
In both cases we discuss the importance of using a damping algorithm. 

© 2012 by the authors. This paper may be reproduced, in its entirety, for non-commercial purposes. 
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1. Introduction 

Hartree-Fock-Bogoliubov (HFB) theory is the method of choice for describing some special 
features of attractive fermionic quantum systems [41] ■ It is a generalization of the famous 
Hartree-Fock (HF) method >35^ used in quantum chemistry. It also generalizes the Bardeen- 
Cooper-Schrieffer (BCS) theory of superconductivity [3] which was invented in 1957 to ex- 
plain the complete loss of resistivity of certain materials at very low temperature. In 1958 
Bogoliubov realized in [5] that the BCS theory was actually very similar to his previous 
works [SI [HI H] on the superfluidity of certain bosonic systems. He adapted it to fermions, 
leading to a model that is now called Hartree-Fock-Bogoliubov, and which is the main in- 
terest of this article. 

In Hartree-Fock-Bogoliubov theory the state of the system is completely determined 
by two operators [3]. The one-particle density matrix 7 is the same as in Hartree-Fock 
theory [34] . whereas the pairing density matrix a describes the Cooper pairing effect. This 
effect can only hold in attractive systems. When the interaction potential is positive, the 
energy is decreased by replacing a by 0. 

One of the most interesting questions in HFB theory is precisely the existence of pairing, 
that is, the non-vanishing of the matrix a for a minimizer. This question has been settled 
in the simpler translation- invariant BCS theory [5l [48j [49l [39 , 18, 23, 26 and in some cases 
for the translation- invariant Hubbard model [3] . But it remains completely open for general 
attractive systems with a few particles, like those encountered in nuclear physics. In |30j . 
the first author of this paper has shown with Lenzmann the existence of HFB minimizers 
for a purely Newtonian system of N fermions, but it is not yet known if a ^ 0. One of the 
purpose of this work is to answer this question numerically. 

The HFB energy is a nonlinear function of 7 and a. Because of nonlinearity, minimizing 
this functional on a computer is not an easy task. The simpler Hartree-Fock model in which 
a = is now well understood from the point of view of numerical analysis [29[ 111) , even 
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though most authors have concentrated their attention to the special case of electrons in 
an atom or a molecule. Cances and Le Bris have studied in [131 112] the simpler fixed point 
algorithm called the Roothaan algorithm [42] and they have shown that this algorithm 
either converges or oscillates between two points, none of them being the solution of the HF 
equation. This result was recently improved by Levitt [31]. Cances and Le Bris have also 
proposed a new algorithm called Optimal Damping, which is now used in several chemistry 
programs. It is based on the fact that one can freely minimize the energy over mixed HF 
states instead of pure HF states, by Lieb's variational principle [33] . 

Because the HFB model is an extension of HF theory, it is natural to believe that these 
ideas can be applied to the HFB case as well. In particular, under appropriate assumptions 
on the interaction potential, Bach, Frohlich and Jonsson have recently shown in [2] an HFB 
equivalent of Lieb's variational principle. Up to some difficulties that will be explained later, 
we will show in this paper that the previously mentioned results can indeed be transposed 
to the Hartree-Fock-Bogoliubov model. 

The paper is organized as follows. In the next section we quickly recall the basic for- 
mulation of Hartree-Fock-Bogoliubov theory. Then, in Section [3j we derive the discretized 
HFB equations and we prove that, in the limit of a large Galerkin basis set, the discretized 
solution converges to the true solution. We also discuss at length the possible symmetries of 
the system and we formulate the theory when these symmetries are taken into account. 

In Section|4]we study the HFB Roothaan algorithm and we prove that it either converges 
to a solution of the HFB equation, or oscillate between two points, none of them being a 
solution of the equations. We then introduce an equivalent of the Optimal Damping Algo- 
rithm of Cances and Le Bris, which is based on an optimization in the set of mixed HFB 
states. 

Section [5] is devoted to the presentation of some numerical results. We first consider 
the purely gravitational model studied by Lenzmann and Lewin [30] and we numerically 
discover that there is always pairing. Then, we introduce a simplified model for nucleons, 
for which we as well present some preliminary numerical results. We particularly discuss the 
importance of using a damped algorithm instead of a simple fixed point method, a fact which 
has already been noticed in nuclear physics [16]. Our approach could help in improving the 
existing numerical techniques. 

Acknowledgement. The authors would like to thank Laurent Bruneau and Julien Sabin 
for useful discussions. They acknowledge financial support from the French Ministry of Re- 
search (ANR-10-BLAN-0101) and from the European Research Council under the European 
Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement MNIQS 
258023). 

2. A quick review of Hartree-Fock-Bogoliubov theory 
2.1. Hartree-Fock-Bogoliubov states and their energy 

We consider a system composed of N identical fermions, described by the many-body Hamil- 
tonian 

JV 

H(N) = J2 T i+ E W *> 

3=1 l^k<e^N 



-1 



Mathieu LEW IN & Severine PAUL 



acting on the fermionic TV-body space Sj N = f\\ f)i where ft is the space for one particle. 
Here, T : S) — > is a one-body operator and W : Sj 2 — > Sy 1 accounts for the interactions 
between the particles. We use the notation Tj for the operator T which acts on the jth 
component of the tensor product Sj N = /\^ $), that is Tj = l(g)---(g)T(g)---®l, and a 
similar convention for Wm- 

Most of what follows is valid in an abstract setting. However, for the sake of simplicity, 
in the whole paper we will restrict ourselves to the special case of nonrelativistic fermions 
with q internal degrees of freedom, moving in M 3 (q — 2 for spin-1/2 particles like electrons). 
We also assume that no external force is applied, and that their interaction is translation- 
invariant. Then, in units where m = 1/2 and h = 1, we have 

i3 = i 2 (K 3 ,C«), T = -A, W ke = W{x k -x e ). 

The TV-body space Sj N = L 2 (R 3 , C q ) consists of wave functions %(xi, o"i, xjy, <7/v) 
which are antisymmetric with respect to exchanges of the variables (a^cr,-). In principle 
W(xk — xg) is also a function of the two internal variables ak,&i G {lj -••> q} of the particles 
k and I. Again for simplicity, we will assume that W only depends on the space variable 
Xk — xg. Finally, we make the assumption that W is smooth and decays fast enough at 
infinity to ensure that H is bounded from below. To make this more explicit, we assume in 
the whole paper that 

W = Wi + W 2 E L P (R 3 ) + L ? (R 3 ) for some 2 < p < q < oo. (2.2) 

Sometimes we will make more precise assumptions on W. 

We are interested in the case where W is attractive iW < 0), or at least partially 
attractive (W $J on a set of measure non zero). By translation invariance, the Hamiltonian 
H(N) has no ground state (that is, the bottom of its spectrum cannot be an eigenvalue). 
But it may have one once the center of mass is removed, if W is sufficiently negative. 

In a nonlinear model approximating the many-body problem above, there could be a 
ground state, even if the system is translation- invariant. Of course, translation invariance 
is not lost and there are then infinitely many ground states, obtained by translating the 
system arbitrarily. In Hartree-Fock theory |35[ |3] , such breaking of symmetry is known to 
occur for instance when W{x) — —\/\x\ is a purely gravitational interaction and T = —A 
(nonrelativistic), or T = y/l — A — 1 (pseudo-relativistic) , see [30j [32] and Theorem 12.31 
below. 

For attractive systems, it is often convenient to allow for another symmetry breaking, 
namely that of particle number. This means that the fixed particle number TV is replaced 
by an operator M whose eigenvalues are 0, 1, 2, .... Only the average particle number is well 
defined for a quantum state. The classical way to define Af is to introduce the fermionic 
Fock space 

which gathers all the possible n-particle subspaces in a direct sum. A (pure) quantum state 
in J 7 is a vector \1/ = ipo ® ipi ® • • • which is normalized in the sense that 

I*£ = l^| 2 + £hM«» = 1 - 
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The average particle number is the diagonal operator 

TV := 0® 0n, 

such that the average number of particles in a state $ is given by the formula 




Instead of imposing that \& £ S) N which is equivalent to *f> being an eigenvector of Af, 
Af^f = N^, we will only fix the average particle number of *f>: 

(*,7V*> = N. 

Allowing to have ip n ^ for n ^ N is useful to describe some physical properties of attractive 
systems. In most practical cases it is expected that the variance J2 n >o( n ~ -^) 2 ll^nlU" wm 
be quite small, i.e. that "J will live in a neighborhood of 9) , 

Similarly to the particle number operator Af, the many-body Hamiltonian H(N) is now 
replaced by a many-body Hamiltonian H on Fock space 

H:=Offi0i?(n) (2.3) 

which is nothing else but the diagonal operator which coincides with H (n) on each n-particle 
subspace. We will not discuss here the problem of defining H as a self-adjoint operator on 
T. 

The Hartree-Fock-Bogoliubov (HFB) model generalizes the well-known Hartree-Fock 
(HF) method and it allows for breaking of particle number in a very simple fashion. The 
method consists in restricting the many-body Hamiltonian H on J to a special class of states 
called Hartree-Fock-Bogoliubov states (or quasi- free states), which are completely character- 
ized by their one-particle density matices [3] . 

Let us recall that a state in Fock space has two one-particle density matrices, instead 
of one in usual HF theory. These are two operators 7 : fj — > S) and a : S) — > fj, which are 
defined by means of creation and annihilation operators by the relations [3] 

(vM^/Ms)*)^ = (.9,7J%, (*,o(/)a(ff)% - (ff, 

When "J lives in a particular iV-particle subspace S) N , then a(f)a(g)^ £ S) N ~ 2 hence 
a(f)a(g)^f) jr — for all /, g £ Sj and the matrix a vanishes. However for a general 
state f e J, one can have a 7^ 0. 

The two operators 7 and a satisfy several constraints. First, we have 7* = 7, ^5 7 ^ 1 
(in the sense of operators) and Tr7 = (^>,Af^>) — N, for the one-particle matrix 7. On the 
other hand, the so-called pairing matrix a satisfies a T = —a. Its kernel a(x, a, x' , a') — 
—a(x',a',x,a) can thus be seen as an antisymmetric two-body wavefunction in Sj 2 . It is 
interpreted as describing pairs of virtual particles, called Cooper pairs. 

It is well known that a quantum state "J £ Sy N such that (7*) 2 = 7* is necessarily a 
Slater determinant (that is, a Hartree-Fock state). The same is true for states in Fock space. 
Consider a pair (7, a) which is such that 

r 2 = r, with r : = ( \ a _) (2.4) 
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on S) © .f). Hence we have for instance aa* = 7 — 7 2 . Then there exists a unique state ^/ in 
which has 7 and a as density matrices. This state has the property that any observable 
can be computed using only 7 and a, by Wick's formula (see Thm 2.3 in [3]). The quantum 
states obtained by considering projections V are called Hartree-Fock-Bogoliubov states and 
they generalize usual Hartree-Fock states. When a = 0, then 7 = Ylj=x IVjXVjl i s a rank- TV 
projection and the corresponding state is the usual Slater determinant 

* = <px A • •• A (fN = ^Jjyj ^et {Pi( x j 1 a j)) ' 

When a^O, >t can be obtained by applying a Bogoliubov rotation to the vacuum but we 
will not explain this further. For the present work, we will only need the formula of the total 
energy, in terms of 7 and a: 

= Tr(-A) 7 +i/ / W{x-y)(p 1 {x)p 1 {y)-\ 1 (x,y)\ 2 + \a(x,y)\ 2 )dxdy 
1 it 3 JR3 v ' 

:=f(7,a) (2.5) 

where p-y(x) — Trc? (7(2;, x)) is the density of particles in the system. The terms in the 
double integral are respectively called the direct, exchange and pairing terms. Taking a = 
one recovers the usual Hartree-Fock energy which has been studied by many authors [35l [38| 
[H [3] . Our main goal in this paper is to investigate the minimization of the more complicated 
nonlinear functional £(7, a), when 7 and a are submitted to the above constraints, and 
its numerical implementation. We will show below that the energy £ is well defined in an 
appropriate function space, under our assumption (|2.2p on W. 

Note that the variance of the particle number for a HFB state \& in Fock space can be 
expressed only in terms of a by 

(*, (M - AAI>)^) 2 *)^ = J> " Nf |^„|||„ = 2 Tr^aa*), 

see Lemma 2.7 in [3]. The spreading of the HFB state over the different spaces Sj n is therefore 
determined by the Hilbert-Schmidt norm of the pairing matrix a. We recover the fact that 
an HFB state has a given particle number if and only if its pairing matrix a vanishes. 



2.2. Pure vs mixed states 

In our previous description, we have only considered pure states in Fock space, that is states 
given by a normalized vector $ £ J. For practical purposes, it is very convenient to extend 
the model to mixed states, which are nothing else but convex combinations of pure states, 
given by a (many-body) density matrix in T 

^ = E A ^->^I with A ^°< E^ = 1 > *<- >l '.-.; 'V 

j j 

The average particle number and energy are then given by the formulas 

Tr^ (MD)=Y1 X i <*i > ) > Tr ^ ( H£) ) = I] X J (*i - »j ) • 
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Resorting to mixed states is mandatory at positive temperature, when the equilibrium state 
of the system will actually always be a mixed state. But it is also very useful at zero 
temperature, even if the true ground state is a pure state. We will recall later the practical 
advantages of using mixed HFB states. 

Similarly to what we have explained in the previous section, there is a class of mixed 
Hartree-Fock-Bogoliubov states, which are completely characterized by their one-particle 
density matrices 7 and a. The latter now satisfy the constraint 

on Sj © Sj, which is nothing else but the relaxation of the constraint T 2 = T of pure states. 
The set of one-particle density matrices of mixed HFB states is therefore a convex set, 
whose extremal points are the density matrices of pure HFB states. One should remember 
that the set of mixed HFB states is not the convex hull of HFB pure states, however. The 
relation between the density matrices (7, a) and the corresponding HFB states in T is highly 
nonlinear. 

The energy of a mixed HFB state described by the density matrices (7, a) is given by 
the same formula (|2.5j) as for pure states. Hence, minimizing this energy under the relaxed 
constraint (|2.6|) is equivalent to minimizing the full quantum energy over all mixed HFB 
states. The natural question arises whether a minimizer, when it exists, is automatically a 
pure state. The answer to this question is positive in many situations, as we will see below. 

Before turning to the comparison between the minimization among pure and mixed 
states, we first introduce the variational sets on which the energy is well defined. The sets 
of all pure and mixed HFB states with finite kinetic energy are respectively given by 

V := {(7, a) € 61(f)) x 62(f)) : a T = -a, T = T* = T 2 , Tr(-A) 7 < 00} (2.7) 

and 

K := {(7, a) e 61 (S) x 6 2 (3) : a T = -a, < T = T* < l fl(M >, Tr(-A) 7 < c»} , (2.8) 

(the matrix T is the one appearing in (|2.4|Q . Here 61(f)) and 62(f)) denote the spaces of 
trace-class and Hilbert-Schmidt operators [33] . The expression Tr(— A)7 is to be understood 
in the sense of quadratic forms, that is 

3 

Tr(-A)7 = ^Tr(p k jp k ) E [0,+oo], with p k = -id Xk . 

k=l 

In practice we want to fix the total average number of particles. For this reason we also 
define the constrained sets 

V(N) := {(7, a) e V : Tr 7 = N} (2.9) 

and 

K(N) :={(j,a) €K : Tr 7 = 7V}, (2.10) 



of pure and mixed states with average particle number N. In practice N is an integer but it 
is convenient to allow any non- negative real number. 
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The following lemma says that the energy is a well-defined functional on the largest of 
the above sets IC, and that it is bounded from below on JC(N) for any TV ^ 0. 

Lemma 2.1 (The HFB energy is bounded-below on K.(N)). When W = W\ + W 2 £ 
L P (M 3 ) + L q (R 3 ) with 2 sc: p ^ q < oo, then £(7, a) is well defined for any (j,a) £ JC. It 
also satisfies a bound of the form 

V( 7 ,a)e/C, £( 7 ,a)>iTr(-A) 7 -C(iV) (2.11) 

for some constant C(N) depending only on N = Tr(7). 

Proof. The assumption that W = W1 + W2 £ L P (R 3 ) + L q (R 3 ) with 2 ^ p ^ q < 00 implies 
that VF is relatively form-bounded with respect to the Laplacian, with relative bound as small 
as we want [15]. This means \ W\ ^ e(— A) + C e in the sense of quadratic forms, for all e > 
and for some constant C c . This can now be used to verify that the energy is well defined 
under the assumption that Tr(— A)7 < 00. First, we have for the direct term 

\W{x-y)\ Pl {x) Pl {y)dxdy^eN f \V^\ 2 + C e N 2 ^ eATr(-A) 7 + C e N 2 , 

JR 3 

where in the last line we have used the Hoffmann-Ostcnhof inequality 27J, 

/ \VVfh\ 2 < Tr(-A) 7 . (2.12) 

JR 3 

The exchange term is bounded similarly by applying the inequality \W\ ^ e(— A) + C e in x 
with y fixed: 

\W(x-y)\ h{x,y)\ 2 dxdy s$ e / / \V xl {x,y)\ 2 dx dy + C t / | 7 (x, y)\ 2 dx dy 

Jr 3 Jr 3 Jr 3 Jr 3 

= eTr(-A) 7 2 + C* e Tr 7 2 ^ e Tr(-A) 7 + C e N, 
since 7 2 ^ 7. Similarly we have, since aa* ^7 — 7 2 ^ 7, 

\W(x-y)\ \a(x,y)\ 2 dxdy < Tr (e(-A) + C £ )aa* < e Tr(-A) 7 + C e N. 



All this shows that all the terms in the energy are well defined when (7,0) £ K.(N). Also, 
we have 

£(7, a) > (1 - e - eN/2) Tr(-A) 7 - C t N - C £ 7V 2 /2. (2.13) 
Taking e = 1/(2 + N) finishes the proof. □ 



Lemma 12.11 allows us to define the minimization problems for pure and mixed states as 
follows: 

UN) := inf £(7, a), (2.14) 

J(N):= inf £(7, a). (2.15) 

Of course we have J(N) ^ I(N) since V{N) C JC{N). In many cases, we have that I(N) = 
J(N) and that any minimizer, when it exists, is automatically a pure HFB state. We give 
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two results in the literature going in this direction. The first deals with purely repulsive 
interactions and it is Lieb's famous variational principle |33) (see also Thm. 2.11 in [3]). 

Theorem 2.1 (Lieb's HF Variational Principle [33j ) . Assume that 

and let N be an integer. Then for any (7, a) € IC(N), there exists (7',0) G V(N) such that 

£(7, a) >£(Y,0). 

In particular, we have I(N) — J(N). 

If W > a.e., then any minimizer for I(N), when it exists, is necessarily of the form 
(Y,0) with {if = 7'. 

We see that for repulsive interactions, W > 0, there is never pairing (a = 0) and the 
ground state is always a pure HF state, that is, a Slater determinant. The fact that, in HF 
theory, one can minimize over mixed states and get the same ground state energy is very 
important from a numerical point of view. This was used by Cances and Le Bris [13| 112] to 
derive well-behaved numerical strategies, to which we will come back later in Section [4.21 

For purely attractive interactions, the following recent result of Bach, Frohlich and Jon- 
sson [2] is relevant. 

Theorem 2.2 (HFB Constrained Variational Principle |2j). Assume that the number 
of spin states is q = 2 (spin-1/2 fermions) , and that W can be decomposed in the form 

W(x- y) = - [ d(i(u))l Aui (x) 1 Alo (y) (2.16) 

on a given measure space with A u a measurable family of bounded domains in R 3 . 

Let N be any positive real number. Then for any (7,0) € K.(N), we have 

£(7,a)>£( 7 ',c/), (2.17) 

with 

7' = 5®(J°), = VsO^ff)® (_ X J) (2-18) 

(the matrices act on the spin variables), and 

_ t _ — _ Trt + tu + 7n + m 

4 

This HFB state is pure: (7', a') € V(N). In particular, we have I(N) = J(N). 

Furthermore, ifW < a.e., then any minimizer for I(N), when it exists, is necessarily 
of the previous form. 

Remark 2.1. Note that N does not have to be an even integer in this result. Since 
TrL 2 (R 3 )(s) = N/2, the operator g must have one eigenvalue different from 1 when N is 
odd, and it follows that a 7^ in this special case. 

In Theorem I2.2[ we have decomposed the operator 7 acting on L 2 (R 3 x {j,4-},C) ac- 
cording to the spin variables as follows: 

hn 74.-A 
\m 7u/ ' 
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Theorem 12 .21 savs that when W satisfies (|2.16j) . one can minimize over states which are pure, 
real, and have a simple spin symmetry. The antisymmetry of a is only contained in the spin 
variables, hence the Cooper pairs are automatically in a singlet state. Of course, one can 
express the total energy only in terms of the real operator g, as follows 

£(j',a') = 2 Tr L2(R3) (-A) 5 

+ f f W{x-y)(2p g {x)p g {y)-\g{x,y)\ 2 + Wg{l-g){x,y)\ 2 ). 

JR3 J R 3 \ / 

In practice it will be more convenient to keep a pairing term a(x,y) not a priori related to 
g and to optimize over both g and a, that is, to consider mixed states. When W satisfies the 
assumptions of the theorem, any ground state will automatically lead to a = ±y/g(l — g). 

Let us conclude our comments on Theorem l2.21 by noticing that several simple attractive 
potentials W can be written in the form (|2.16p . For instance the Fefferman-de la Llave 
formula [T7] 

1 1 f°° dr f 

shows that a simple Newtonian interaction W(x — y) = —\x — is covered (here ^B(z,r) 
is the characteristic function of the ball centered at z, of radius r). Hainzl and Seiringer 
showed in [25) that any smooth enough radial function W can be written in the form 

W{x-y) = I drg{r) I dz t B lz,r)(x) 1 B(z,r) (y) 
Jo Jm? 

for some explicit function g, whose sign can easily be studied. 



2.3. Existence results and properties of minimizers 

Before turning to the discretization and the numerical study of the HFB minimization 
problem, we make some comments on the existence of minimizers. In addition to the infinite 
dimension and the nonlinearity of the model, an important difficulty is the invariance under 
translations. For instance, there are always minimizing sequences which do not converge 
strongly (assuming there exists a minimizer, one can simply translate it far away). 

Consider the electrons in an atom or in a molecule, with fixed classical nuclei (Born- 
Oppenheimer approximation). From the point of view of the electrons, the problem is 
no more translation-invariant, once the nuclei have been given a fixed position. Since the 
Coulomb interaction between the electrons is repulsive, W(x — y) = \x — y\~ x , Theorem 12. II 
tells us that there is never pairing, a = 0. In this case there are several existence results, 
starting with the fundamental works of Lieb and Simon [35] and continuing with works by 
Lions [38], Bach [I], Solovej [46l [47]. 

For interactions W which have no particular sign, the pure HF problem was studied by 
Friesecke in [19] and by the first author of this paper in [32] . There is an HVZ-type theorem 
for HF wavefunctions which states that some binding conditions imply the existence of 
minimizers (Theorem 22 in [32]). 

After the fundamental paper of Bach, Lieb and Solovej [3] with its study of the Hubbard 
model, to our knowledge the existence of ground states for the HFB model with pairing 
was only studied recently by Lenzmann and the first author of this paper in [30 . Some 
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caricatures of HFB in nuclear physics had been previously considered by Gogny and Lions 
in [22]. 

We give here an existence result for the variational problem I{N) . In some cases we have 
I(N) = J(N) (see the previous section) and for this reason, we only concentrate on I(N). 
The next theorem can be proved by following the method of |30j , which dealt with the more 
complicated case of a pseudo-relativistic kinetic operator yl — A — 1. 

Theorem 2.3 (Existence of minimizers and compactness of minimizing se- 
quences). We assume as before that W = Wx + W 2 £ L p {R 3 )+L q (R 3 ) with 2 p ^ q < oo. 
Let A > 0. Then the following assertions are equivalent: 

(1) All the minimizing sequences (7„,a„) C IC(X) for J(A) are precompact up to translations, 
that is there exists a sequence (xt) C R 3 and (7, a) € IC(X) such that, for a subsequence, 



lim 



(l-A)^(r afc7nfc 7-_ XJk - 7 )(l-A) 



- AW 2 



lim 

k— > 00 



61 

(l-A)V 2 (r Xk a nt 



In particular (7, a) is a minimizer for /(A). 
(2) The binding inequalities 

I(X) < I(X - M ) + I(fi) 

are satisfied. 



for all < fi < X 



= 0. 



(2.19) 



Furthermore, if W is Newtonian at infinity, that is 

W(a0<-r7 for a > and \x\ ^ R, 
\x\ 



(2.20) 



then the previous two equivalent conditions are verified. 

The assumption that the interaction is Newtonian at infinity is a big simplification, as 
it means that two subsystems receiding from each other always attract at large distances. 
One can expect that minimizers exist even if the potential is not attractive at infinity, as 
soon as it has a sufficiently large negative component. A typical effective potential W(x) 
used in nuclear physics is nonnegative for small |x|, and has a negative well at intermediate 
distances |41j . At infinity it typically decays like for two protons, and exponentially 

fast when one of the two particles is a neutron. Even in HF theory, we are not aware of 
any existence result dealing with such potentials, however. We will numerically investigate 
a model of this form in Section 15.31 

The form of the nonlinear equation solved by minimizers is well-known in the physics 
literature, and it was re-explained in [3]. The following result summarizes some known prop- 
erties. 

Theorem 2.4 (HFB equation and properties of minimizers [3, 30j). A HFB mini- 
mizer on JC(N) solves the nonlinear equation 



1 



(-00,0) 



[F r -hN)+S 



where ^ 8 1{ }(-Pr — A^V) h as ^ e same form as T, and where 

1 \ A 7 
v -1 ' r U* -h- 



(2.21) 



(2.22) 
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with h~ ( = —A + p 7 * W — W(x — y)~f(x, y) and n(x, y) — a(x, y)W(x — y). 

IfW(x — y) = —k\x — J/l^ 1 (Newtonian interaction) and N is an integer, then all the 
minimizers are of the special form (|2.18[) . In this case, we have either a = and 7 is a 
projector of rank N , or a 7^ and 7 has an infinite rank. 

The nonlinear equation (|2.2ip is in principle similar to the usual equation obtained in 
HF theory, 

7 = l(-co, M )(M+<5 (2-23) 

Indeed, (|2.21[) reduces to (|2.23|) when a = 0. Let us however emphasize that the mean-held 
operator has a spectrum which is symmetric with respect to 0. Hence Fr is usually not 
even semi-bounded, on the contrary to h 7 which is always bounded from below. Furthermore, 
the operator Af does not commute with Ft (except when a = 0) and the equation cannot 
be written in a simple form as in HF theory. This will cause several difficulties to which we 
will come back at length later. 

The fact that 7 has an infinite rank when there is pairing, a ^ 0, is a dramatic change 
of behavior compared to the simple HF case. However, no information on the decay of the 
eigenvalues of 7 seems to be known. 

An important open question is to show that minimizers actually exhibit non-vanishing 
pairing a ^ 0, at least for a sufficiently strong attractive potential W. On heuristic grounds, 
one expects such a phenomenon of "Cooper pair formation" to be energetically favorable due 
to the attractive interaction among particles. However, it seems to be a formidable task to 
find mathematical proof for this claim. The existence of pairing is known in some particular 
situations (when N is odd and W is Newtonian, see Remark [2~Tl for the Hubbard model [3], 
or in translation-invariant BCS theory [5lli8 l l39 l l39 lH8l[23l l26]). but for the model presented 
here, we are not aware of any result of this sort. One of the purposes of this paper is to 
investigate this question numerically. 

3. Discretized Hartree-Fock-Bogoliubov theory 

In this section, we write and study the Hartree-Fock-Bogoliubov model in a discrete basis. 



3.1. Convergence analysis 

Here we show that the ground state HFB energy in a finite basis converges to the true 
HFB ground state energy when the size of the basis grows. We consider a sequence of finite- 
dimensional spaces Vn C H 1 (R 3 , C 9 ) for h ^ 0. We assume that any function / £ H 1 (R 3 , C q ) 
can be approximated by functions from Vh'- 

V/ G H\R 3 , C«), 3f h G V h such that ||/ - f h \\ m — > 0. (3.1) 

h— fO 

We typically think of a sequence Vh given by the Finite Elements Method. Let iTh denote 
the orthogonal projection on Vh in L 2 (R 3 ,C q ). We define the set of density matrices living 
on Vh (with average particle number N) as follows 

fch(N) = {(7, a) G JC(N) : 717,7717, = 7, 717, a Wi = a). (3.2) 

The corresponding minimization problem is now 

I h {N)= inf 5(7, a). (3.3) 
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Since K.h{N) C K(N) by definition, it is obvious that Ih(N) ^ I(N). The following result is 
a consequence of Theorem 12.31 



Theorem 3.1 (Convergence of the approximate HFB problem). When W — W\ + 
W2 G L P (R 3 ) + L 9 (R 3 ) with 2 $5 p ^ q < 00 and under Assumption (|3.1| on the sequence 
(Vh), we have 

lim Ih(N) = I(N). (3.4) 

h— s-0 

// i/ie binding inequality (|2. 19|) is satisfied, then any sequence of minimizers (7^,0/1) G 
Kh{N) for Ih{N) converges, up to a subsequence and up to a translation, to a minimizer 
(7, a) € K-(N) of I(N), in the sense that 



lim 



(l-A) 1 / 2 (r Xt7 / lfc ^- 7 )(l-A) 1 / 2 



lim 



©1 

(l-A)^ 2 (T Xk a hk r^ Xk -a) _ =0. (3.5) 

e 2 



Proof. We only have to show that Ih(N) —> I(N) as h — > 0. Then, any sequence of exact 
minimizers (7^, ah) for Ih{N) is also a minimizing sequence for I(N). Applying Theorem l2.3l 
concludes the proof. 

We know that finite-rank operators are dense in K.(N). Let (j,a) € IC(N) be any such 
finite rank operator. Let (fi)fLi be an orthonormal basis of the range of 7. By Lowdin's 
theorem (Lemma 13 in |32j). we know that the two-body wavefunction a can be expanded 
in the same basis (/1, fx)- Now, for every i — 1, .., K, we apply (|3.1|i and take a sequence 
ft G Vft, be such that f£ — > fi in H 1 (M. 3 ). The system is not necessarily orthonormal 

but we have — > (fi,fj) = Sij- Applying the Gram-Schmidt procedure, we can 
therefore replace the {f^fLi by an orthonormal set (gi)fLi C 14 having the same properties. 
An equivalent procedure is to take g\ = ^2^—i{S h 2 )jif^ where Sh is the Gram matrix 
{(fi 1 fj 1 ))- Let now Uh be any unitary operator on L 2 (R 3 ) which is such that Uhfi = 9i for 
all i = 1, K. We then take jh := UhjU^ and a/j := UhaU^ . In words, we just replace the 
fi by in the decomposition of 7 and a. To see that (7^, ah) G K,(N), we just notice that 

lh a h \ = fU h _0_\ / 7 a \ /[/,: _0 

Ofci-w v° V«* 1-7/ v 

Note also that Tr(~fh) = Tr( 7 ) = N since Uh is unitary. Now 7^ and belong to K.h(N) by 
definition, hence we have that £(jh,&h) Ih(N). On the other hand, by the convergence 
of f' 1 (hence of <?- 1 ) towards fi in iJ 1 (R 3 ), we easily see that 

lim£ ('jhjCth) = £(l> a ) 

by continuity of £. Therefore we have proved that 

lim sup l h (N) ^ £ (7, a). 

/i-s-0 

This is valid for all finite-rank (7, a) G K{N), hence we deduce that 
lim sup /ft (TV) < inf £(7, a) = I(N). 

h^O (7,«)6?C(JV) 

On the other hand the inequality Ih(N) > -f(^V) is always satisfied, hence we have proved 
the claimed convergence Ih(N) — > I(N). □ 
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3.2. Discretization 

In this section we compute the HFB energy £(7, a) of a discretized state (7, a) G fCh(N) 
and we write the corresponding self-consistent equation. We fix once and for all the ap- 
proximation space Vh and we consider a basis set ixi)f=i 01 Vk> which is not necessarily 
orthonormal. We will assume that Vh is stable under complex conjugation, which means 
that / € Vh => / € Vh (this amounts to replacing Vh by Span(V/ l , Vh)). Then we can choose 
a basis (Xi)i=i which is real, that is Xi — Xi f° r au * — lj— J-N&. This will dramatically 
simplify the calculation below. 

Since i^hl^h — 7 and TThaWh = a, we can write the kernels of 7 and a as follows 

N b JVi, 

l[ x i V)a,a' = X! G vXi( x )<yXj{y)a', a(x,y) ay(T , = ^2 A ijXi{x)aXj(y)(r'- (3-6) 

The complex conjugation on %y is superfluous but we keep it to emphasize the difference 
between 7 and a. The matrices G and ^4 (defined by the previous relation) satisfy the 
constraints G* — G and A T = —A. Note that since A is antisymmetric, we can also write 

<x(x, y)„ >a > = A a {xi{x)aXj{y)v - Xj(x)aXi{y)v') 

= V2 Y A tjXt^Xj(x,cr,y,(j') 

l<i<j<JV(, 

where (xi A Xj)(x, o~, y, a 1 ) := ( 

Xi{ x )aXj (y)<r> - Xi(y)a'Xj( x )a)/V2 is a two-body Slater 
determinant. Let us also remark that (|3.6[) can be written in the operator form 

JV,, N b 

7 = Gl 3 \X-i){Xj\, " = Y Al 3 IXi)<XJ|- 

Again the complex conjugation on \~j is superfluous. 

Let us define the so-called overlap matrix S = S* = S by 

1 r « r 

Sij = (Xi,Xj)sj = z2 Xi(x)aXj(x)<7 dx = V / x l (a;) -x i (a;) .da;. (3.7) 

It is tedious but straightforward to verify that the constraint 



,0 QJ \a* 1 - 7/ \Q 1 

can be written for the matrices G and A in the form 



0\ /SGS SAS \ /S N 

M 0,1*0 c c7=?c J 1 n c J 



v ooy"V<sA*s s-sgsj " Vo 

or, equivalently, 

sC TST < T (3.9) 

where T and S are defined by 
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Another way to write the constraint is ^ S 1 ^ 2 TS 1 ^ 2 ^ 1. Let us notice that we have used 
everywhere the fact that S = S* = S = S T . The formulas are much more complicated when 
S is not real. 

The energy can be expressed in terms of the matrices G and A, as well. A calculation 
shows that 

a) = Tr(hG) + ~ Tr(G J{G)) - \ Tr(GK(G)) + \ Tr(A* X(A)). (3.11) 

The trace here is the usual one for Nb X Nb matrices. As we think that there is no possible 
confusion with £(7, a), we will also denote by £ (G, A) this discretized energy functional. In 
formula ([3~TT]) . 

9 f 



N b N b N b 

J(Gh = £ {ij\tk)G m K{G) l3 = (ityj)Gkt, X(A) tj = ^ (ik\jt)A u , 

k,l=l k,£=l k,£=l 

(3.12) 

and 

(ij\k£):= f f W(x-y) X i(xy Xj (x) X k(yT Xi(y)dxdy. (3.13) 

We use here the notation a*b = X^=i ^ b a (but the complex conjugation is superfluous for 
our real basis, hence K = X). Similarly, we have 

Tr( 7 ) = Tr(SG). (3.14) 

We define the discretized number operator as 

N=(?°) (3.15, 



v S, 

The constraint Tr( 7 ) — N can be written equivalently as 

Tr(NT) = 2N — Nb 

We deduce from this calculation that the variational problem Ih(N) can be written in 
finite dimension as 

I h (N) = min \S{G, A) : TST < T, Tr(NT) = 2N - 7V 6 | (3.16) 

where we recall that T and S have been defined in f|3 . 10|) . Here the infimum is always 
attained because the problem is finite dimensional. 

In this form, the discretized problem is very similar to the usual discretized Hartree-Fock 
problem [HIES], in dimension 2Nb instead of Nb. There is a big difference, however. In HF 
theory the constraint involves the matrix S instead of N. This difference will cause several 
difficulties. To understand the problem, let us introduce a new variable T' = S 1 ^ 2 TS 1/ ' 2 
(which is the same as orthonormalize the basis (Xi))- Then the constraint ^ TST ^ T is 
transformed into T' ^ 1. However, the constraint on the number of particles becomes 
Tr(S" 1/2 NS" 1/2 T') = 2N - N b . In usual Hartree-Fock theory, the matrix S~ 1/2 NS _1/2 is 
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replaced by the identity. The fact that this matrix then commutes with the Fock Hamiltonian 
(defined below) simplifies dramatically the self-consistent equations. Here we get 

s -l/2 NS -l/2 = 

which commutes with the Fock Hamiltonian if and only if A = 0. 

The self-consistent equation is obtained like in [3J. The result is as follows. 

Lemma 3.1 (Discretized HFB equation). Let T be a minimizer for the variational 
problem Ih(N). Then there exists fj, € K such that T solves the linear problem 

min | Tr(F T - /iN)f : s$ f ST «C T j (3.17) 

where 

Fx := X ^ , h G := h + J(G) K(G). (3.18) 

The solution T can be written in the form 

T = S~ 1/2 Ut-oo.o) (S~ 1/2 (F T - ,uN)S- 1/2 ) + S^j S~ 1/2 (3.19) 

where ^ S ^ 1 lives only in the kernel o/S _1 ^ 2 (Fx — /iN)S -1 ^ 2 . 

The solution T of the self-consistent equation f|3 . 19[) may be equivalently written by 
considering the generalized eigenvalue problem 

(F T - ,iN)/ 4 = a S f u (f h Sfj) = Sij. (3.20) 

Then we have simply (assuming ^ for all i = 1, 2Nb) 

£,<0 

Again, this is similar to the Hartree-Fock solution [TT] [55] except that /i is unknown and N 
does not always commute with Fx- 

Remark that although the basis functions Xi are reai ) the density matrix T is not neces- 
sarily real. In the next section, we will restrict ourselves to real- valued density matrices and 
impose some spin symmetry. 

3.3. Using symmetries 

The HFB energy £(T) has some natural symmetry invariances which we describe in detail 
in this section. Recall that since £ is a nonlinear functional, it cannot be guaranteed that 
the HFB minimizers will all have the same symmetries as the HFB energy. The set of all 
minimizers will be invariant under the action of the symmetry group, but each minimizer 
alone does not have to be invariant. 

We have already allowed for the breaking of particle-number symmetry and we hope to 
find an HFB ground state. It will then automatically break the translational invariance of the 
system. There are three other symmetries (spin, complex conjugation and rotations) which 
are of interest to us. We have the choice of imposing these symmetries by adding appropriate 
constraints, or not. Because this reduces the computational cost, it will be convenient to 
impose them. 



Co-i) 
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3.3.1. Time-reversal symmetry 

Let us now assume that q — 2, which means that our fermions are spin-1/2 particles. Since 
the Laplacian and the interaction function W do not act on the spin variable, the HFB 
energy has some spin symmetry, which can be written for q = 2 as 

Vfc = i,2,3, £(E fc rs*) = £(r) 

where 

'io k 



St- := 



ia k 



with (Ti, 02,03 being the usual Pauli matrices. Note that £fc has the form of a Bogoliubov 
transformation, hence EfcT££ is also an HFB state. The number operator is also invariant 
which means that 

£ fe A^ = Af 

for all k = 1, 2, 3. Thus, the contraint Tr(7) = N is conserved and we have E^TE^ € /C(JV) 
when r € /C(AT). 

Another important symmetry is that of complex conjugation which means this time that 

£{T) = S(T) 

and which is based on the fact that the Laplacian and W are real operators. Again we have 
Tr(7) = Tr(7) hence K.{N) is invariant under complex conjugation. 

As was explained in [24] (see Remark 5 page 1032), the density matrices (7, a) can be 
written in the special form 

7' = g ® ( p ^ , a' = a <g> ( ^ ^ , with g = g T = g and a = a T = a (3.21) 

(the 2x2 matrix refers to the spin variables) , if and only if 

|E fc rE^=r, for A =1,2, and (g ^ 

In other words, T is invariant under the action of the group generated by Sx, E2 and C. This 
invariance is sometimes called the time-reversal symmetry. As remarked in |24| . imposing 
SfeLE^, = F for all k = 1,2,3 implies a = which is not interesting for us. 

When W can be written in the form (|2.16[) . the Theorem 12.21 of Bach, Frohlich and 
Jonsson [2], tells us that there is no breaking of the time-reversal symmetry. That is, we can 
always minimize over such special states. For other interactions W this is not necessarily 
true but it is often convenient to impose this symmetry anyhow. 

Because it holds 

-F^rs* = Efci^rSfc, F T = F r , 

it can then be verified that minimizcrs under the additional symmetry constraint, satisfy 
the same self-consistent equation as when no constraint is imposed. 

When we discretize the problem by imposing time-reversal symmetry, we use two real 
symmetric matrices G and A, related through the constraint that 

< TST < T, with T := \ and S = ( S Q ° V (3.23) 
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The energy becomes 

£(G, A) = 2Tr(hG) + 2 Tr(G J(G)) - Tr(G K (G)) + Tr(AK(A)). (3.24) 

and the associated particle number constraint is Tr(SG) = N/2. In practice we always as- 
sume that N is even for simplicity. The basis (xi) is now composed of (real- valued) functions 
in i? 1 (R 3 ,IR), instead of functions in H 1 (M. 3 ,C q ) as before, and the formulas for S, h, J, K 
and X are the same as before. 

3.3.2. Rotational symmetry 

The group SO (3) of rotations in K 3 also acts on HFB states and it leaves the energy invariant 
when the interaction W is a radial function. In this section we always assume that the spin 
variable has already been removed according to the previous section and we denote by g 
and a the corresponding (real symmetric) density matrices. If the spin were still present, 
rotations would act on it as well. 

To any rotation R € SO (3) we can associate a unitary operator on L 2 (R 3 ), denoted 
also by R, defined by (Rip)(x) = (p(R~ 1 x). Now we say that an HFB state V with density 
matrices (7, a) is invariant under rotations when it satisfies 

&T@* =T, where ^ = 

Note that & is a Bogoliubov rotation since R = R. The density matrices of an invariant 
state satisfy 

g(Rx,Ry) = g(x,y), a{Rx, Ry) = a(x,y) 

for all x, y G K 3 and any rotation R e 50(3). 

As the angular momentum L = x x (— iV) generates the group of rotations, a (smooth 
enough) HFB state is invariant under rotations if and only if 

J^r = rj^, where 5? = 
The density matrices g and a can then be written in the special form 

= ■^^29 e (\x\,\y\)(2e+l)Pi(u x -Lj y ), a{x,y) = -L ^ a e (\x\, \y\) (2£+l) P t (u x -u y ) 

(3.25) 

where Pi is the Legendre polynomial of degree £, which is such that P((l) — 1. The constraint 
on g and a are transfered in each angular momentum sector (labelled by t ^ 0), leading to 

(00\ < (9 t a* W10\ 

on L 2 ([0, 00), r 2 dr)©L 2 ([0, 00), r 2 dr). However, there is no such constraint between different 
^'s. The total average particle number is given by 

Tr(. 9 )=^(2£+l)Tr(/)=7V/2. 
This is the only constraint which mixes the different angular momentum density matrices. 
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Now we can discretize the radial HFB problem. We choose a finite-dimensional subspace 
VJ-ad in i 2 ([0, oo), r 2 dr) with basis (%i, XN b ), which we use to expand the density matrices 
g and a . Then we fix a maximal angular momentum £ max and we truncate the series 
in (|3.25l) . This is the same as taking as discretization space 



V = 



{f(\x\)Y e m (u x ) : / e Kad, £ max , -l^m^l] C L 2 (K 3 ,M) 



where Y^ m is the spherical harmonics of total angular momentum £ and azimuthal angular 
momentum m. We then assume that g and a are radial and live in this space. The matrices 
G l and A 1 of and a f in the basis (xi) are defined similarly as before by 

N b N b 

g%r') = £ Gf i Xi(r)x J -(r'), O = £ 4 x^V)- (3.26) 

»,i=i »,j'=i 

The constraints on the matrices G^ and ^ are 

G £ A e \ , /SO 



with 



< T ST < V, with T' := ( A , _ G A and S=l (3.27) 



Sij = / Xi(r)Xj{r)r 2 dr 
Jo 



and 



J2(2£+l)Tr(SG e ) = N/2. (3.28) 

£=0 

The total HFB energy is now 
S (G°, G^ max , A , A'-") =2£(2U1) Tr(//G £ ) 

^ (2^+ 1)(2/ + l)(2Tr(G £ J(G £ ')) - Tr(G £ i^'(G £ ')) + Tr(A* , (3.29) 



+ 

e,t'=o 
where 



and 



Kj ■= / X'iir) Xj (r) r 2 dr + £(£ + 1) / Xi (r) X j (r) dr, 
■/o Jo 

iV b JVi, 

J(G*)ij:= £ (y|nm) ,oG< n , K*(G?) ir .= g (im\jn)i^> G l mi 

rn , n — m.n— 



/>OG /'OC 

{ij\mn)i t i> '.= I r 2 dr s 2 ds Xi(r)Xj(r)Xm(s)Xn(s)w e ^(r,s) 
Jo Jo 



w t ,t> (r, S )-=\J W (Vr 2 + s 2 - 2rst) P t (t) P t > (t) dt. 
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Any minimizer (G°, G fmax , A , A 4 ""**) of £ under the constraints ([3~27]l and (|3~28| will 
be of the form 

T £ = ]T /f (/f) T , < * < Cax, (3.30) 

ef<0 

where the vectors /^'s solve the generalized eigenvalue problem 

(F e - fi(2e + 1)N) ft =e*S ft, (f u Sfi) = 5 l3 (3.31) 

with 

h t q\ £ — (2J(G 1 ')-K U \G 1 ') K U '(A 1 ') \ 

-h l ) + ' -2j(G r ) + X"'(G £ ') J ' 

The Euler-Lagrange multipler /i appearing in (|3.31|) is common to all the different angular 
momentum sectors and it is chosen to ensure the validity of the constraint ()3.28j) . 

4. Algorithmic strategies and convergence analysis 

In this section we study the convergence of two algorithms which can be used in practice 
to solve the HFB minimization problem (|2.15[) . In order to simplify our presentation, we 
restrict ourselves to the finite-dimensional case, that is, to the discretized problem ()3.16j) . 
We also assume that the discretization basis (xj) is orthonormal, such that S = hN b , the 
(2A;,) x (2Nb) identity matrix. Finally, we only consider states which are invariant under 
time-reversal symmetry like in Section 15. 3. II This means that the HFB state is described by 
real and symmetric matrices G and A such that 

°°)^T:=( G A (4.1) 

ooy " \a 1-gJ ^ v° v 

The energy is given by (13.241) . 

£{T) = 2Tr(/iG) + 2Tr(G J(G)) -Tr(GK(G)) + Tr(AK(A)). (4.2) 

The extension to more general situations is straightforward. 

The energy £ is continuous (it is indeed real- analytic) with respect to T. Also the set 
K, of density matrices T of the form (|4.ip is compact in finite dimension. Hence minimizers 
always exist and, as we have seen, they solve the nonlinear equation 



l(-oo 



, 0) (F T -/iN) +5, (4.3) 



where /i is a Lagrange multiplier chosen to ensure the constraint that Tr(G) = N/2. Of 
course we must have N/2 ^ iV&, the dimension of the (no-spin) discretization space Vh, 
otherwise the minimization set is always empty. 

4.1. Roothaan Algorithm 

The most natural technique used in practice to solve the equation (|4.3|) is a simple fixed 
point method [40|, 116] . This is usually refered to as the Roothaan algorithm in the chemistry 
literature [42]. The iteration scheme is the following 



T„ +1 = l ( _ 00 ,o) F T „ - Mn+iN +<W (4-4) 



A Numerical Perspective on Hartree-Fock-Bogoliubov Theory 



21 



At each step, one has to determine fi n +i such as to satisfy the constraint Tr(G„+i) = N/2. 
If the operator Fy„ — /i rt +iN has a trivial kernel, then 5 n+ i = 0. This is the usual situation 
encountered in practice. In the iteration (I4.4p . the state is assumed to be pure at each step 
(T is an orthogonal projection). Recall that by Theorem 12.21 we know that minimizers of 
£ under a particle number constraint are always pure, under suitable assumptions on the 
interaction potential W . The algorithm is stopped when the commutator 

l|[ T «^T„]|| 

and/or when the variation of the HFB state 

are smaller than a prescribed e. 

Our purpose in this section is to study the behavior of the Roothaan algorithm (14. 4p . 
In the Hartree-Fock case, it was shown in a fundamental work of Cances and Le Bris [T2l 
113] . that the algorithm converges or oscillate between two points, none of them being a 
solution to the equation (|4.3[) . This result was recently improved by Levitt in [31]. We will 
explain that the results of Cances-Le Bris and Levitt can be generalized to the HFB model. 
Actually, in a discretization space of dimension Nf,, HFB is equivalent to a Hartree-Fock- like 
minimization problem in dimension 2iVf,, with additional constraints. The adaptation of the 
previously cited works in the HF case reduces to handling these contraints. 

In order to avoid the convergence problems of the Roothaan algorithm, Cances and Lc 
Bris have proposed the Optimal Damping Algorithm (ODA). We will study the equivalent 
of this algorithm in HFB theory in Section 14.21 

To start with, we show that the Roothaan algorithm is well defined, in the sense that for 
any HFB state T„, there exists (T n+ i, (J, n +i, S n+ i) solving (|4.4p . To this end, we follow [12l 
113] and introduce the auxiliary functional 

£(T, T') Tr(hG) + Ti(hG') + 2 Tr(G J{G')) - Tr(G K{G')) + Tr(A K{A')) (4.5) 

as well as the variational problem 

7 T (A) := min {f(T, T') : Tr(G') = a} (4.6) 

which consists in minimizing over T' with T fixed. The matrix T' must be an admissible 
HFB state which, in our context, means that 

(o o) « ' - {% i -g) < G ' (G ' )T - V - G ' {A ' )T = W = A '' (47) 

Recall that we have chosen an orthonormal basis and that the spin has been eliminated. It is 
clear that (|4.6j) admits at least one solution T', as soon as ^ A ^ JV&, where we recall that 
Nf, is the dimension of the discretization space Vh- The following says that these solutions 
are exactly those solving the equation of the Roothaan method. 

Lemma 4.1 (The Roothaan algorithm is well defined). The function A € [0, iVf,] i-» 
Jx(A) is convex, hence left and right differentiable. For any A € [0, the minimizers T' 
of -fx (A) are exactly the states of the form 



T' = 1(-oc,o)(Ft-m'N)+<5' 



(4.8) 
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where // € [d-Iy(X),d + I r (X)] and 0^5'^ 1 {0 }(F T - fi'TSS). If $ er(F T - n'N), then 
§' = and T' is unique, for any such (if G [9_Jx(A), &fJx(A)]. 

Proof. To see that lx(A) is convex, let ^ Ai < A2 ^ Nj, and let be a minimizer for 
J T (Ai) with i = 1,2. Then iT^ + (1 -i)T!, is a test state for the problem I T (iAi + (1 -t)X 2 ). 
Therefore it holds 

It (i Ai + (1 - t)X 2 ) ^ £(T, tT[ + (1 - t)T' 2 ) = t£(T, T[) + (1 - t)£(T, T' 2 ) 

= tI r (X 1 ) + (l-t)I r (X 2 ). 

Then, by convexity we get that It (A') > /x(A) + /j(A' — A) for any A' € [0,-/V&] and any 
(i e [d-I r (X),d + I r (X)\. Thus 

I T (A) - fiX = min{/ T (A') - /jA' : sC A' < 7V h } 
= min{£(T,T')-MTr(G')} 

= Tr(hG) + - minTr(F T - uN)T'. 

2 T' 

In the previous two mins, T' is varied over all possible HFB states, without any particle 
number constraint. It is well known that the minimizers of the problem on the right side are 
exactly the solutions of the equation (|4.8I) . □ 

Lemma 14.11 tells us that for any given T„, there always exists at least one solution 

(T 

n+ij fJ-n+i j <5n+i ) of the equation (|4.4|) . It is obtained by solving the minimization problem 
Iy n (N/2), and one has to take ^ n +i <= [d-Ir n (N/2),d+Ir n (N/2)]. We can always take by 
convention 

_ cL/ Tn (N/2) + d+I Tn (N/2) 
Hn+i ■— 2 ' 

However, T„ +1 is not uniquely defined yet because of the possibility of having 5 n+ i 7^ 0. As 
we have seen it is unique when 

0^(F T „-/j ll+1 N). (4.9) 

In this section we always assume that it is possible to find (T„ + i, /-i n +i, 5n+i) exactly. In 
practice, we will only know (T„ + i, S n+ i) approximately. Later in Section l4~3l we explain 

how to do this numerically. We will also see that the condition (|4.9[) is "very often" satisfied. 
This vague statement is made precise in Lemma T4.4I below. Following Cances and Le Bris, 
we now introduce the concept of uniform well-posedness. 

Definition 4.1 (Uniform well-posedness). We say that for a given initial HFB state 
T , the sequence (T n ) generated by the Roothaan algorithm is uniformly well posed when 
there exists 77 > such that 

|F T „ - Mn+iN| > n (4.10) 

for all n > 0, where fi n+1 = (<9_ J Tn (iV/2) + d+I Tn (N/2))/2. 

Note that the condition |Fx„ — /i n +iN| ^ 77 is equivalent to (— 77, 77)ner(Fx„ — // n +iN) = 0. 
Later in Section [4.31 we will make several comments concerning Assumption (|4.10p . 
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We have seen that the sequence generated by the Roothaan algorithm can be obtained 
by solving the minimization problem 

Iy n = min£(T„, T'). 
T' 

Since £ (T, T') = £(T' , T), we conclude that the Roothaan algorithm is the same as mini- 
mizing £ with respect to its first and second variables one after another, inductively. This 
fact allows to prove the following result, which is the HFB equivalent of Theorem 7 in [T5] 
and Theorem 5.1 in [31] in the HF case. 

Theorem 4.1 (Convergence of the Roothaan algorithm). Assume that < N/2 < 
Nb. Let To be an initial HFB state such that the sequence (T ra ) generated by the Roothaan 
algorithm is uniformly well posed. Then 

• The sequence £ (T2„, T2 n +i) decreases towards a critical value of £ ; 

• The sequence (T2„,T2„+i) converges towards a critical point (T,T') of£; 

• IfT = T' ; then this state is a solution of the original HFB equation (|4.3[) . but ifT ^= T', 
then none of these two states is a solution to (|4.3[) . 

Theorem 14.11 says that (provided it is uniformly well posed) the sequence T„ will either 
converge to a solution of the self-consistent Equation (|4.3[) , or oscillate between two points 
T and T', none of them being a solution to the desired equation. 



Proof. We split the proof into several steps. 



Step 1: /i„ is uniformly bounded. It will be very useful to know that the sequence 
u n is uniformly bounded. The following says that, as soon as we fix Tr(G) = N/2 with 
< N/2 < Nb, the chemical potential /i cannot be too negative and too positive. 

Lemma 4.2 (Bounds on the multiplier /i). Let T' be any fixed HFB state and 



T M :=l ( _ OOi0 )(F T ,-|iN) 



G> A^ 
A^ 1-G A 



the corresponding HFB ground state at chemical potential /i. There exists a constant C which 
is independent of T' and a, such that 

VH^-C, TrG.^^r and V> ^ C, Tr G M ^ iV b - - . (4.11) 

The lemma says that the average number of particles in l(_oo.o) (Fx — A*N) tends to Nb 
when fi — > oo whereas it tends to when /i — > — oo, this uniformly with respect to the state 
T' used to build the mean- field operator Fx' • 



Proof. We first remark that there exists a constant G such that 

IFt'KC (4.12) 

for any HFB state T'. This follows from the fact that F~c< is continuous with respect to T' 
and that the latter lives in a compact set since we always have ^ T ^ 1. The chosen norm 
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for I-Ft'I does not matter since we are in finite dimension. Now, for fj, large enough we can 
use regular perturbation theory and obtain that 



|l(-oo,o)(Fr' -A*N) -l ( -oo,o)(-/iN)| 



(-00,0) 



Ft' h 



Note that 



l(-co,o) I -rrN| = < 



Taking the trace against N gives the result. 



10" 

K 0, 

'o N 
,0 1, 



AN I - l(_oo,0) I -rrN 



for /x > 0, 
for [i < 0. 



C_ 

R 



□ 



From Lemma 14.21 we deduce that our sequence /i„ is bounded. Indeed, since we have by 
construction Tr(G„+i) = N/2 with < N/2 < Nb, we must have 



max C 



2C_ 
~~N 



^ Hn+i ^ max C 



C 



N b - N/2 



(4.13) 



otherwise Tr(G„ + i) would be too small or too large. 



Step 2: convergence of £(T n , T rl +i). We now follow jT2l H31 [10] . At each step, we know 
from Lemma |4 . 1 1 that T n+1 is a solution of the minimization problem miny £(T n , T'). In 
particular, we deduce that 



£ (T n , T„ + i) ^ £ (T„, T„_i) — £ (T„_i, T ra ). 



(4.14) 



Thus the sequence of real numbers £ (T„, T n+ i) is non-increasing. It is also bounded from 
below, hence it converges to a limit t. We now use the uniform well-posedness to prove an 
inequality which is more precise than (|4. 14[) . We remark that 

£ (T„, T„_i) — £ (T„, T n+ i) — — Tr Fx„ (T n _i — T n +i) 

= 1 Tr (F Tn - Mn+1 N) (T w -i - T n +i) 
^ i Tr |F T „ - Mn+ iN| (T„_i - T„ +1 ) 2 
^ I Tr (T n _x - T n+1 ) 2 = I IT^x - T n+1 || 2 . (4.15) 

In the above calculation we have used that Tr N T„ + i = Tr N Y n _i = N — Nf,. We have also 
used that T„ + i is the negative spectral projector of Fx„ — jtz n +iN, such that we can write 



Mn+1 N ) 



■ /X„+lN (T n+1 - T„ + i) 



Finally, we have used that 7 < 1 is equivalent to (7 - P) 2 < P- 1 (7 - P)P ± - P(j - P)P, 
for any orthogonal projector P, thus 

n+iy^n-l — 1 n+1 J J- „+l ~ J- n+1 { J- n-1 ~ i-n+l) i-n+1 \± n-1 — 1 n+1 ) ■ 



A Numerical Perspective on Hartree-Fock-Bogoliubov Theory 



25 



Since £ (T ra , T n+ i) converges to a limit £, we deduce that 

y"l ||T n+ i - T n _i| 2 < oo. 

In particular, ||T n+ i — T n _i|| — > which is called numerical convergence in [12l [13] 



Step 3: convergence of T2 n and T2n+i • In order to upgrade the numerical convergence 
to true convergence, we use a recent method of Levitt [31]. Namely we show that 

(f(T„,T n _i)-*) 9 - (£(T n ,T„ +1 )-£) 8 ^ITn-i-Tn+il (4.16) 

for a well chosen < 9 ^ 1/2. Summing over n and using the convergence of £(T n , T„_i), 
hence of {£ (T n , T„_i) — then gives the convergence of T 2n and T 2li +i. 

For the proof of (|4.16[) . we argue as follow. Consider a (real, no-spin) pure HFB state T. 
It is possible to parametrize the manifold of pure HFB states around T by using Bogoliubov 
transformations as follows: 

H ^ e H Te- H 

where H is assumed to be of the form 

h P\ uT T . t — 

, h = -h = -h, p =p = p. 

-p -hj 

These constraints ensure that iH is a self-adjoint Hamiltonian such that e H = e~ %l S" H ' is a 
Bogoliubov rotation. They also ensure that e H Te~ H stays real. That H i-> e H Ye~ H is a 
local chart of the manifold of pure HFB states around T follows from the arguments in [3] 
as well as simple considerations in linear algebra. 

Let us now consider the energy £ in a neighborhood of any fixed (T, T'). The map 

/ : (H,H') ^ £(e H r e - H ,e H 'r'e- H ') - ^ Tr Ne H Te~ H - | TrNe ff ' T' e~ H ' 

is real analytic in a neighborhood of (0, 0) for any fixed and any fixed pure HFB 

states (T, T'). The Lojasiewicz inequality (Theorem 2.1 in [31]) then tells us that there exist 
< 9 ^ 1/2 and a constant k > such that ||iJ|| + ||iJ'|| ^ n implies 

\f(H,H') - fm 1 - 9 < k-^IVh/^^OI + \Vwf(H,H')\ 
A simple computation shows that 
V H f(H, H>) = ~ [F r , - M 'N , e ff Te- ff ] , V fl - f(H, H') = I [F T - , e^re^'] . 

If we rephrase all this in our setting, this means that for any fixed pure HFB states (Ti, T^) 
and any fx, fj.' € R, there is a k > such that for any (T 2 , T 2 ) another pure HFB state which 
is at most at a distance k from (Ti, T^), we have 

£{T 1 ,r[) -f(T 2 ,T 2 ) + M 'TrN(G 2 -G 1 ) + ^TrN(G 2 -G' 1 )| 1 ~ e 

< k- 1 ( || [F T , - M 'N , T 2 ] || + || [F T2 - M N , T' a ] || ) . (4.17) 
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The constants k and 6 depend on /x, fi 1 and of the reference point (T^T^). But they stay 
positive as soon as /i, fjf and (Ti,T^) stay in a compact set. By a simple compactness 
argument, we therefore deduce that there exists a neighborhood of the compact set 

|(T,T') : £(T, T') = £, Tr(G) = Tr(G') = iV/2 j 

such that for any (T, T') in this neighborhood and /x, /x' in a compact set in R, we have 

1-0 

£(T, T') -g + n'(N/2- TrNG) + /i(iV/2 - Tr NG') 

^ k- 1 ( || [F T , - //N , T] || + || [F T - , T'] || ) (4.18) 

for some < 9 ^ 1/2 and some k > 0. We recall that I is by definition the limit of 
£ (Y n , T„ + i) . 

Recall our inequality (|4. 13|) which says that /i„ is uniformly bounded. Also, we know that 
£(T n , T„ + i) converges to £ so, for n large enough, (T„, T n+ i) must be in the neighborhood 
of the level set I. Choosing \i = (i n +i and \j! = \i n and using that G„ and G n +i have the 
correct trace, we get the estimate 

(S (T„, T n -i) - i) ^ < n- 1 (|| [F Tn - Mn+1 N , T tt _i] || + || [PT n _, - MnN , T„] ||) 

= kT 1 || [F T „ - m«+iN , T n _i - T„ + i] || 
^ G ||T n _i — T n+ i|| 

for n large enough. Here we have used that T„ commutes with Ft ii _ 1 -/i„N and that T Jl+ i 
commutes with Fx„ -/J n +iN by construction, and that ||Fx„|| and /i n +i are both uniformly 
bounded. In order to conclude, we use the concavity of x i-» x 6 and (|4.15|) like in [31] to 
obtain 



(£(T„,T n _!) - (£(T n ,T„ +1 ) 



t 



^ — ^ x _q T n _i) — £(T„,Y n+1 ) 



2(£(T„,T„_ 1 ) -£) 

>rfl/{2C) |T n+1 - T^rl 

by (|4.15|) . This concludes the proof of the inequality (|4.16[) . hence the proof of the conver- 
gence of (T2n, T2n+i), towards some pure HFB states (T, T'). 

Step 4: the limit (T, T') of (T2n,T2„+i) is a critical point of £. Since we have 
T~2n — > T and T 2n +i — > T', we deduce that Fx 2 „ — > Fx and Fx 2n+1 — > Fx', by continuity of 
the map T i-> Fx - Extracting a subsequence, we can assume that H2n k - > A*' an d H2n k +i 
We have 

Y2n fc = l(-oo,0)( F T 2 „ fc _ 1 - M2n fc N), T 2 „ fc + i = l(-oo,0) ( F Y 2nfe - M2n fe +1 N ) 

and, by uniform well-posedness, 

| F T 2ni .-i - M2n fc N| ^ T], |F T2 „ t - Ai2« fc + lN| ^ 7]. 
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Passing to the limit k — > oo we get 

T = l ( _ OO!0) (F T ,- M 'N) and T' = l ( _ oo , 0) (F T - fiN). 

This exactly means that (T, T') is a critical point of £ on Vh(N /2) x Vh(N/2). Note that 
we have also |Fx' — m'N] ^ 77 and Fx — A«N| ^ 

The remaining statements are verified exactly like in the HF case. This concludes the 
proof of Theorem 14.11 □ 



4.2. Optimal Damping Algorithm 

In the previous section we have studied the convergence properties of the Roothaan algo- 
rithm, which consists in solving the self-consistent equation by a fixed point method. We 
have seen that the algorithm can either converge or oscillate between two states, none of 
them being a solution to the problem. 

Examples of such oscillations in quantum chemistry have been exhibited by Cances and 
Le Bris |T2l [13] . In this case the potential W is repulsive and there is no pairing. In order 
to cure this problem of oscillations, Cances and Le Bris proposed in 13J a relaxed algorithm 
called the Optimal Damping Algorithm (ODA). This method makes use of the important fact 
that one can minimize over mixed states and get the same ground state as when minimizing 
over pure states only (Theorem 12. II) . 

The same oscillations can a priori happen in HFB with an attractive potential W . 
They are frequently seen with the Roothaan algorithm and we will give several numerical 
examples later in Section[5] Even when the sequence T„ eventually converges towards a single 
state T, these oscillations can slow down the convergence considerably. This phenomenon 
is well known in nuclear physics. Decharge and Gogny already advocate in [16] the use of a 
damping parameter between two successive iterations, in order to "slow down the convergence 
on the density matrix. In this way the average field varies slowly and we can insure the 
convergence on the pairing tensor step by step" (see [16j page 1574). Even in the modern 
computations, this damping parameter is fixed all along the algorithm (Nathalie Pillet, 
private communication) . 

We suggest to transpose the method of Cances and Le Bris to the HFB setting by using 
an optimal damping parameter, chosen such as to minimize the energy. This means resorting 
to mixed states even if the final ground state is always a pure HFB state. This is theoretically 
justified when the assumptions of the Bach-Frohlich-Jonsson Theorem 12.21 are fulfilled. 

The ODA involves two density matrices T„ and T„. The HFB state T„ is always pure 
but T n can (and will usually) be a mixed HFB state. The starting point To = To being 
chosen, the sequence is then constructed by induction as follows: 

(1) One finds (T n+ i, fJ, n +i) solving 

T n+1 = 1(^00,0) (Ff n ~ Mn+iN) and Tr(G„ +1 ) = A/2. 
This is always possible, by Lemma 14. II and we can take as before 

d-I fn (A/2) + d+I fn (A/2) 

Mn+l : = " 2 ' ' 

in case is in the spectrum of Fx„ — /i„+iN. 

(2) One lets 

T„+i = t n+ iT n + (1 — t„+i)T„+i 
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where the damping parameter t n+ i € [0, 1] is chosen such as to minimize the (quadratic) 
function 

t^£(tf n + (l-t)T n+1 ). 

(3) The algorithm is stopped when || [T„, Fx„] || and/or ||T n+ i — T„| are smaller than a 
prescribed e. 




iN) 



Fig. 1. The Optimal Damping Algorithm of Cances & Le Bris in the HFB case. 

The general strategy of the ODA is displayed in Figure [T] By construction we see that 
£ (T„) is a non-increasing sequence. This guarantees the convergence of the ODA. The result 
is the following 

Theorem 4.2 (Convergence of the ODA). Assume that < N/2 < N b . Let T = f 

be an initial HFB state such that the sequence (T„) generated by the ODA is uniformly well 
posed, that is 

Vn, |F tn - Mn+iN| ^ r] > 0. (4.19) 

Then 

• The sequence £(T n ) decreases towards a critical value of £ ; 

• The sequence T„ numerically converges towards a critical point T of £ , in the sense 
that T„ + i — T„ — > 0, T n+ i — T„ — > and that all the limit points T of subsequences of 
(T„) solve T = l ( _ oo ^ 0) (F T - /iN). 

Proof. The proof is exactly the same as in the Hartree-Fock case [lUJ [14] and we only 
sketch it. First we have by definition £(T„ + i) ^ £(T„), so £ (T„) must converge to a limit 
I. Now we have 

£(T„+i) = £((1 — iri+i)T rl + t„+iT n+ i) = £{T n ) — t n+ ia n+ i + tf l+1 b n+ i 

where 

t n+ i = argmin te[0 x] ( - ia„+i + £ 2 £> rl+ i) 
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and with 



a n+ i := Tr (T„ - T„+i) = Tr |F^ - /xN| (f „ - T n+i ) ^ ?/ 



b n +i — 2Tr(G„+i — G n )J(G n +i — G n ) — Tr(G„+i — G n )-ftT(G n +i — G n ) 

+ Tr(i n+1 - A n )K(A n+1 - A n ) 



In finite dimension we have |b n +i| ^ C 
prove that 



Y n +i — T„ 



^ (C/rj)a n+1 . This can be used to 



— tn+ia.n+i + t n+1 b n +i ^ — ea n +i ^ — en 
for some e > independent of n. This now proves that 



T rl — T ra +i 



Tl+1 



< oo, 



hence that T„ — T„ + i — » 0. In order to conclude the proof, we notice that 
T n +i — T n = T n +i — T n + (1 — t n ){in-i — T n j 

which finally implies 

II 2 

^ ||T„ - T„ + i|| 2 < oo and ^ T„ - T n+ i < oo. 

n n 

Since T, i+ i = l(~oo.o)(Ff — Mn+iN) by definition, the proof that any limit T of a subse- 
quence of (T n ) satisfies the self-consistent equation is elementary. □ 



4.3. Handling constraints 

Both the Roothaan algorithm and the ODA are based on Lemma 14.11 which says that for 
any given Fx, there exist //, 8' and T' such that 



T' = l ( _ OOi0) (F T -//N)+ ( 5', 
TrNT' = N ~N b . 



(4.20) 



The purpose of this section is to explain how to solve this problem numerically. To simplify 
our notation, we consider in this section a generic matrix 



h p 
P ~h ) ' 



with p = p = p and h = h = h 



and we study the problem consisting in finding T, /i and 6 such that 

fT = l ( _ oo , )(F-MN)+5, 
[TtNT = N - N b . 

Assume first that p = (Hartree-Fock case) . Then we have 

h 
-h 



(4.21) 



(4.22) 



30 



Mathieu LEW IN & Severine PAUL 



which commutes with N. The solution of (|4.22p is then given by the aufbau principle, 

1- G )' = 1^^ + 6 

where [i is the (iV/2)th eigenvalue of h, counted with multiplicity and S lives in the corre- 
sponding eigenspace. Equivalently, 

K K' 

G = Y^ v i v f + Ul ViV i 

i=l i=K+l 

where the v^s solve the eigenvalue equation 

hvi = ti Vi, 

K = Tr l(-oo,ejir/2) C 1 ) i s * ne dimension of the direct sum of all the eigenspaces corresponding 
to the eigenvalues < epj/2 an d K' = Trl^^ >£N , 3 ](h) is the dimension of the direct sum of 
all the eigenspaces corresponding to the eigenvalues ^ £n/2- The »j's are chosen such that 

£ N 

o s$ m < i, k+ 22 n * = y • 

Therefore, finding T, /i and 5 in the Hartree-Fock case only requires to diagonalize h once. 

In the Hartree-Fock-Bogoliubov case (p ^ 0), the situation is more complicated since N 
does not commute with F. Let us consider the real function 

^ < \ TrNl ( _^ 0) (F- M N)+iV b 
W ■ V- ^ v{n) = . (4.23) 

We are interested in solving the equation 

up(fx) = N/2. 

In the Hartree-Fock case, Vp is a non-decreasing piecewise constant function. There is a 
solution /i to vf(h') — N/2 when N/2 belong to the range of vp. Otherwise, one has to 
partially fill a shell using the matrix S. 

In the Hartree-Fock-Bogoliubov case, vp is also non-decreasing and in general it is much 
smoother when p ^ 0. The following lemma summarizes some important properties of vp in 
both the HF and HFB cases. 

Lemma 4.3 (Elementary properties of v). Let F be as in (|4.21|) . Then the function 
Vp defined in (|4.23j) is increasing with respect to \i. It can only have finitely many jumps. It 
satisfies for some constant C depending only on 

• v(fi) ^C/n for [i ^ — C; 

• up(fi) ^N b - C/fi for fi^C. 

I/O i o-(F-/iN), then 

% ) = 2 VMM^ (4.24) 

e 3 >0 

where (F — /iN)«j = EjUj. 

Proof. The behavior of vp for 3> 1 was already studied in Lemma T4. 2 1 
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The matrix F — fiN is a linear function of fi € R, hence by [28], we know that its 
eigenvalues form a set of real analytic functions. They cannot be constant because the 
matrix N does not vanish. The eigenvalues of F — /iN all behave like ±fi for large fi, by 
perturbation theory. We conclude that can be an eigenvalue of F — fiN for a finite number 
of //'s, say /Zi < • •• < /Uif. On the other hand, the map fi H ► l(_oo,o) (F — ^N) is real-analytic 
outside of the /i^'s (see [28 J. So j/p is itself real-analytic outside of this set and it can have 
at most a finite number of jumps. 

Outside of the fik's, it is possible to compute the derivative of i/p by usual perturbation 
methods [28] . The answer is (|4.24[) and the fact that dvp jdfi ^ proves that vp is increasing 
with respect to fi, in between these points. That the jumps are all positive can be easily 
seen by an approximation argument using Lemma 14.41 below. We skip the details. □ 

The shape of the function vp is very different in the HF and HFB cases. For a Hartree- 
Fock state, the function vp is piecewise constant and it has jumps at the eigenvalues t\ < 
■ ■ ■ < £N b of ha- The size of the jumps is equal to the multiplicity of the associated eigenvalue. 
An HFB state will most always have a very smooth vp. Of course, the smaller p in the 
Hamiltonian F, the more vp looks like a step function. 

In Figure [2] below, we show the function vp for different values of the pairing term. More 
precisely, we have randomly chosen two symmetric real matrices h and p of size Nb = 5, and 
we display the function vp when the pairing is replaced by tp for t = (Hartree-Fock case), 
t = 0.1 and t — 1. Figure [3] is a plot of the eigenvalues of F — /iN for t = 0.1, as functions 
of fi. Note that there are some crossings of eigenvalues above and below the real line (recall 
that the spectrum is symmetric with respect to 0). But, around the crossings are avoided 
and there is a gap. 

If we repeat the numerical experiment with several random matrices h and p, we never 
see any jump for vp. The purpose of the next result is to clarify this observation. 

Lemma 4.4 (Generic behavior of vp). The Fock matrixF — fiN is invertible if and only 
if h±ip — fi are invertible. More precisely, 

miner(|F - fiN\) = min M|(/i + ip - m) _1 || _1 , \\(h - ip - /i) -1 ^ 1 ) . (4.25) 

The set of real symmetric matrices h and p such that 

cr(F - /LtN) n {0} = for all fie M. 

is open and dense in {(h,p) : h = h T = h, p = p T — p} . For h and p in this set, vp is 
real-analytic on R. 

It is obvious that there are matrices h and p for which F — fiN has as eigenvalue for 
some fi e R. The simplest examples are HF Hamiltonians for which p = and F — fiN is not 
invertible each time fi equals an eigenvalue of h. If p does not vanish but commutes with h, 
then we have \h + ip — fi\ 2 = \h — ip — fi\ 2 = (h — fi) 2 + p 2 and we see that is never in the 
spectrum of F — fiN when the kernel of p does not contain the eigenvectors of h. However 
there are counterexamples with p invertible not commuting with h. For instance, F + N is 
not invertible for 




We now turn to the proof of Lemma l4~4l 
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Fig. 2. The function vy{p) which gives the average number of particles in the state 1(— oo,0)(F ~ I'N), in 
terms of the chemical potential fi. The pairing term in F is equal to tp with t = (Hartree-Fock case, red 
curve), t = 0.1 (blue curve) and t = 1 (green curve). 




-10 12 3 

Fig. 3. The eigenvalues of F — in terms of fj. for t = 0.1. 



Proof. The operator F — /iN is unitarily equivalent to 

(ii) (F -" N »(4?)^u--/ + »^) 



From this we deduce that F — /iN is invertible if and only if h + ip — fi and h — ip — jj, are 
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invertible. Then we have 

(F-//N)- X = 
and 

" (F - AiN 




-l 



= max 



( || {h + ip-p) 1 \\ , \\(h-ip-p) 



The statement now follows from the fact that, on a dense open set, the spectra of h±ip do 
not intersect the real axis. □ 



Lemma 14.41 is interesting when we apply the Roothaan or the ODA, because it means 
that, as soon as p ^ 0, most often we will have no choice for /U n +i and we will take <5„+i = 0. 
Saying differently, it is really reasonable to assume that the sequence generated by the 
Roothaan and the ODA are uniformly well posed (of course when the final state is believed 
to have a non vanishing pairing), as we did in Theorem 14.11 and 14.21 

Even if it is in general smooth, the function v can still vary quickly and this will be 
the case when the pairing term p is small. The appropriate method to find the solution of 
^(/f) = N/2 then depends on the properties of i^f- If the Hamiltonian F has a large enough 
pairing matrix p, then i/f is smooth and we can use a Newton-like method to solve the 
equation isp{p) — N/2. A trial chemical potential pP being given, we compute the derivative 
dv-p/dp{p°) using Formula (|4.24[) and then let 



^:= M V(A72-Mm°))(^°)) 



The method can be iterated until convergence of p n towards the desired fi. The convergence 
is very fast, as soon as v-p is smooth. 

If the Hamiltonian F has a small pairing matrix p, the function j/p will be smooth but 
close to a step function. Its derivative varies very quickly and the previous Newton method 
is not appropriate. In this case we can use a simple bisection method. The bounds on vf(p) 
for large \p\ can be used to find a good starting interval [pi, p r ] such that up{p{) < N/2 and 
M^r) > N/2. 

We have to find a new p n +i a t each step of the Roothaan or ODA. It is of course not 
efficient to find /U n +i with a very high precision all along the algorithm. Decharge and Gogny 
advice in Section II. E of [16] to apply the Newton scheme only once at each step. This means 
that 



where v n is the function v corresponding to Fx = Fy„ . This is then the same as doing 
perturbation theory on first order. We use a slightly different strategy which we explain in 
the next section. 



5. Numerical results 

In this section, we present some numerical results for two very simple interactions W . We 
start by considering in Section 15.21 a purely 3-dimensional gravitational model in which 

W(x) = -± g>0. 
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Then we consider in Section 15.31 a (repulsive) Coulomb potential which is perturbed at 
intermediate distances by an attractive effective potential, as is usually employed in nuclear 
physics: 

W(x) = -j— | — a\ exp (-6i|x| 2 ) + a 2 exp [~b 2 \x\ 2 ) , n > 0. 
\x\ 

In the next section we quickly explain our numerical technique to treat these two simple 
systems. 

5.1. Method 

To simulate our physical systems, we have used the open source software Scilab [33]. Our 
potential W is always real, radial and spin-independent. To reduce the numerical cost we 
have therefore always imposed the spin, time-reversal and spherical symmetry. This means 
that we have to cope with £ max + 1 real and symmetric Nf, x Nj, matrices G and A (the 
one-particle density matrix and the pairing density matrix in the ^th angular momentum 
sector). The total energy of the system is given by Equation (I3.29|) and we have to impose 
the constraints fl3.27[) and (|3.28p which we recall here for convenience: 

< T^T* ^ T\ with T e := (J g _f_ Q )j and S = g °) , (5.1) 

!„,„ 

J2(2e+l) r Tr(SG e ) = N/2. (5.2) 

We choose a simple basis set (xi> —iXNt) of L 2 ([0, oo), r 2 dr), made of "hat functions" 
associated with a chosen grid 

= r < 7*1 < • • • < r Nb < r Nb+ i := r max . 

We impose Dirichlet boundary conditions at r max . We have tested different types of grids 
and there was no important difference between them. The results presented here are all with 
regular grids. As we will explain later, for a given basis size Nf,, the results usually depend 
a lot on the value of the radius r max of the ball in which the system is placed. 

Our main goal is to investigate the existence of pairing. We therefore always start by 
doing a precise Hartree-Fock calculation, for which we use the Optimal Damping Algorithm 
described in Section 14.21 We take as initial state a simple uniform state 

N 

G init = 2 Trfffj IdjVb 

and we run HF until convergence. We have observed a global stability of the results with 
respect to initial states, hence the previous simple choice is appropriate (but more clever 
choices might decrease the total number of iterations) . Then, we use the converged HF state 
G op t as initial datum for the HFB algorithm. Of course we have to perturb it a little bit 
since any HF solution is also an HFB solution. We proceed as follows. Assuming that the 
overlap matrix S = Id^ and that £ max — for simplicity, the optimal HF state G op t can 
be written in the form 

N/2 

Gopt = J2 VkV * ' 
fc=l 
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where v k are the N/2 first eigenvectors of the mean-field matrix h, 

hv k =e k v k - 

We then choose a number n v of valence orbitals and a mixing parameter 8, and we perturb 
G op t as follows 

N/2-n v N/2 N/2+n v 

GLit= v kvI + v k vl + {\-6) VkV k> 

fe=l k=N/2-n v + l k=N/2+l 

N/2+n v 
k=N/2-n v + l 

In most cases, we have observed that n v = 1 and 8 = 0.95 works perfectly well, that is, the 
algorithm escapes from the HF solution G opt and converges towards an optimal HFB state. 
But other values of n v and 9 seem to work fine also. 

When the maximum angular momentum £ max is larger than 0, we often first run the 
algorithm with £ max = for a few iterations before switching to the actual value of < m „. We 
stop the algorithm when the commutators [F%, T^] are smaller than a prescribed error. We 
know from (|3.30p and ()3.31j) that these commutators must all vanish for an exact solution 
of the discretized HFB minimization problem. In terms of the matrix S, the right quantity 
to look at is 

X;||s^(FfXS-STfX)S^| 

£=0 

where ||-|| is the usual operator norm for (27Vf,) x (2Nt,) matrices. There is a similar formula 
in the HF case [IT] . 

As we have explained, in the HFB case, ensuring the constraint (|5.2[) is not as easy as 
in the HF case. In the beginning of the algorithm, our state T is rather close to an HF 
state by construction. Therefore, the function vf{h) defined in Section B~B1 is close to a step 
function. We choose an error e and look for the next states T^ +1 having a total number of 

particles Y^i^q(^ + 1) ^(^^n+i) c l° se to N/2, within the error e, using a simple bisection 
method. We use the bisection for a fixed number of global iterations. Then, when the pairing 
term is large enough, we switch to a Newton method in order to find the state T„ + i. We 
have observed that even if in the beginning several Newton iterations can be employed 
at each step, usually only one Newton iteration is necessary after a while. To guarantee 
a good value of the average number of particles in the end, we decrease the error e on 
I Efco x ( 2 ^ + !) TT (SG e n+1 ) - N/2\ along the algorithm. 

5.2. Pure Newtonian interaction 
5.2.1. Model 

Here we consider a system of N spin-1/2 neutral particles, only interacting through the 
Newtonian interaction 

W{x) = —fL, g>0. (5.4) 
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This potential is strongly attractive at short distances. Since l/\x\ does not decay too fast 
at infinity, it is also quite attractive at large distances. The kinetic energy does not scale the 
same as the potential energy. By a simple scaling argument, we can therefore always assume 
that 

9=1- 

This model can be used to describe neutron stars and white dwarfs when N 3> 1. It 
has been particularly studied from a theoretical point of view in the pseudo-relativistic case 
where the kinetic energy is given by T — \J c 4 m 2 — c 2 A — mc 2 , see [361 1371 130] . In our 
simulations we restrict ourselves to the non-relativistic case of the Laplacian T — — A /(2m) 
(in units such that m — 1/2). It would be interesting to take N large but this is of course 
much too difficult from a numerical point of view. 

As mentioned before, we always impose the spin and time-reversal symmetries, which is 
perfectly justified for the ground state since the interaction (|5.4j) satisfies the assumption of 
the Bach-Frohlich- Jonsson Theorem 12.21 We also impose spherical symmetry which, on the 
contrary, is not known to hold for the true ground state. 

One advantage of the Newtonian interaction (|5.4j) is that the operators J and K u can 
be explicitely computed in the basis of hat functions. We have shown in Section f3.3.2l that 
the energy can be expressed in terms of 

/•OO 

„2 , / 2 



(ij\mn) L i, = I r dr I s ds \i(r) Xji r ) Xm(s) Xn(s) w e ,e>(r, s) 
Jq Jo 

where 

m,e(r, s) = \ [ W Ur 2 + s 2 - 2rst) P e (t) P t (t) dt = -\ [ .^f'^ dt. 

A J-i v ' £ j_ x y/ r ^ -f s ^ — 2rst 

Using the well-known formula 

1 min(r, s) 7 



Vr 2 + s 2 - 2rst ^ max(r, s)™+ 1 " W 

v n— v ' 



we deduce that 



lA/f 1 \ min(r, s) n 

2 ^o^- 1 / max(r, s)" +1 

The integral over the Legendre polynomials is related to the usual Clebsch-Gordan coeffi- 
cients as follows 

lf_ i p n (t)w)Pe(t)dt=( £ 1; 2 

and only a finite number of terms are non zero in the sum over n. The final result can be 
expressed as 

(ij\mn)i t £> 



00 / n nl \ roc poo ■ / \ 

i 1 n \ f „2 , / j. mm(r,s) 



n=0 



Oj J rdr J 8 dS max(r,^+i »(r)x*(r)Xm(-)x«(-). (5-5) 



These integrals can be explicitely computed for hat functions and ^ £,£' ^ Imax with £ mayi 
not too large. In our numerical experiments we have put the explicit formulas in Scilab for 
f ra = 1. The integrals were stored in memory during the whole calculation. 
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5.2.2. Roothaan vs ODA 

In the HF case, we have observed that the Roothaan algorithm very often oscillates between 
two states, none of them being the solution of the problem (as described in Theorem 14.11 and 
in [13]). The Roothaan algorithm seems more well behaved in the HFB case. With the model 
presented in this section, we never got real oscillations for HFB. Sometimes the convergence 
is improved by using the ODA, but in most cases the Roothaan algorithm always converges 
towards the same state as the ODA in the end. As we will see later, the situation is very 
different for the model studied in Section [5. 3[ which is inspired of nuclear physics. 

We start by comparing Roothaan and ODA in the HF case. There, oscillations seem 
to be related to the size of the gap between the largest filled eigenvalue and the smallest 
unfilled one. Indeed, oscillations in HF seem to only occur when there is pairing in HFB, 
an effect which is also well-known to be related to the size of the gap (see, e.g., Theorem 5 
in [2]). When there is no pairing, the HF Roothaan algorithm always behaves like the ODA. 
However, the situation is complex and there is no exact rule. Sometimes the Roothaan 
algorithm does not oscillate even when the gap is rather small and there is pairing. 

In Figure U] we display the value of the energy obtained along the algorithm for the 
Roothaan and the ODA, for the following choice of parameters: N = 6, Nf, = 200, £ max = 
and r max = 30. The ODA converges in about 17 iterations, whereas the Roothaan algorithm 
oscillates. We also show the value of the norms \G n — G n -\\ and \G n — G71-2I along the 
Roothaan algorithm. The oscillation between two points is clearly demonstrated. 

When we decrease the parameter r max but keep iVj, = 200 constant, the gap is seen to 
increase slightly and the Roothaan algorithm behaves better. In Table [1] we give the numer- 
ical value of the last filled eigenvalue and the corresponding gap. The Roothaan algorithm 
slowly converges for r max = 25 and it coincides with the ODA when r max = 20. The gap for 
'"max = 20 is 2.5 times the one for r max = 30. We will discuss the occurence of pairing in 
terms of the parameter r max in the next section. 



^max 


£ Af/2 


e N/2+l ~ e N/2 


behavior of HF Roothaan 


20 


-0.532430 


0.159430 


fast convergence 


25 


-0.536706 


0.081016 


slow convergence 


30 


-0.529200 
-0.548554 


0.061928 
0.067422 


oscillations 



Table 1. Value of the last filled eigenvalue an d the corresponding gap ejv/2+1 — e JV/2 m HF, for N = 6, 

Nj, = 200 and ^ ma x = 0. For r ma x = 30 the Roothaan algorithm oscillates and we display the last filled 
eigenvalue and the gap for the two states. 



As we have mentioned the Roothaan algorithm is usually much more well behaved in the 
HFB case. However, sometimes the convergence can be improved dramatically by using the 
ODA. In Figure [5] we display the energy along the iterations of the algorithm in both the 
Roothaan and ODA cases, for Nb = 500, N = 16, £ max = 1 and r max = 10. In this case the 
Roothaan algorithm is very badly behaved. It passes very close to the HF ground state and 
it takes it a very long time to escape from it. On the other hand, the ODA does not suffer 
from this problem and it converges much more rapidly. 
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Fig. 5. HFB energy along the iterations of the Roothaan (blue) and the ODA (red) for JVj, = 500, N = 16, 
^max = 1 and r m ax = 10. The optimal HF energy is also displayed (green). 
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5.2.3. Numerical evidence of pairing 

Pairing effects in a finite discretization basis might depend on the properties of the basis. 
As we have expiained, the occurence of pairing is related to the size of the gap in HF theory 
and this gap varies with the radius r max in which the system is confined. If r max is decreased 
the system is more condensed and the HF gap increases. 

In Figured] we display the HF and HFB ground state energies computed for N = 16 in 
a basis set of size Nb — 200, in terms of r max . The HF and HFB curves are distinct for r max 
large enough and they merge at r max = 6 approximately. This observation is confirmed by 
the value of the norm of A plotted on the right of the same figure. The minima of the HF and 
HFB ground state energies are attained at about r max ~ 10 in the HF case and r max ~ 10.5 
in the HFB case, which is sufficiently far from the merging point. The minima of these curves 
correspond to the best possible approximation for a given basis size Nb (here Nb = 200) and 
a given type of grid (here regular). The difference between the corresponding HF and HFB 
energies is significant. The HF ground state energy at r max = 10 is —19.232176 (in our units 
in which m = 1/2 and e = 1), whereas the HFB ground state energy at r max = 10.5 is 
— 19.240176. The norm of the pairing matrix A is rather large at this point: 

\\A\\ = y/Tr(SAoSAo) + 3 Tr(SAiSAi) ~ 0.462129. 

This goes in favour of the conclusion that pairing really occurs for N = 16 in this model. 
This intuition is confirmed by a more precise calculations with 2V& = 500 which we discuss 
below. 

The observation of pairing requires to have an appropriate r max but it does not require 
to have a very large basis set. Even for Nb = 30 and r max = 10.5, we already find that the 
HFB energy is approximately —19.078416 whereas the HF energy is about —19.072954. The 
corresponding norm of the pairing matrix A is \\A\\ ~ 0.424124. 

Pairing is a subtle effect which decreases the energy by a small amount (much less than 
one percent here) . Catching this effect requires to be very careful when choosing the radius 
r m „. Taking r max too small might lead to the conclusion that there is no pairing. In our 
simulations we have always observed the occurence of pairing, but provided we choose r max 
appropriately. The values of r max at which the HF and HFB energies attain their minimum 
were always found on the right of the merging point of the two curves. In Table [2] below we 
give our results for Nb — 200 and N = 6, 10, 16 and 20. The HFB ground state energy is 
always smaller than the HF energy. 

In the paper [30] , Lenzmann and Lewin have rigorously studied the gravitational model 
of this section. They showed the existence of a ground state in both the HF and HFB cases. 
But, so far, no proof that pairing occurs has been provided. The numerical results of this 
section tend to show that there is actually always pairing, at least for N not too large. 

5.2.4. Properties of the HFB ground state 

In Table[2] below we give our results for Nb — 200 and N = 6, 10, 16 and 20, for the optimal 
values of r max . With £ maK = 1 we have observed that the shells are filled alternatively. In 
HF theory, the cases N = 10 and N = 16 correspond to closed shells, whereas for N = 6 
and N = 20 the last shell is only partially filled. This is a simple explanation for the fact 
that the pairing matrix is much bigger in these cases. 

In Table [3] we display the occupation numbers for the optimal HFB ground state in the 
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7 8 9 10 11 12 13 




Fig. 6. Value of the ground state HF and HFB energies (top) and of the norm of the pairing matrix A 
(bottom), as functions of the radius r max in which the system is confined, for N = 16, iVj, = 200 and 

^max — !• 
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closed shell case N = 16 and in the open shell case N — 20. Because of the spin, these are 
the eigenvalues of Go multiplied by 2 and that of G\ multiplied by 6. Even in the closed shell 
case N = 16, a rather important pairing effect is observed between the last filled orbital 
(the second £ = 1 eigenvalue) and the first unfilled one (the third I = eigenvalue). This 
results in a decrease of the last occupation number of the HF one-particle density matrix by 
approximately 0.228. 



N 


^max 


HF gap 


HF energy 


HFB energy 


\\A\\ 


6 


15 





-1.7327688 


-1.9934252 


1.0242134 


10 


11 


1.023642 


-6.7911634 


-6.8148576 


0.5871951 


16 


10 


1.404396 


-19.232177 


-19.2403096 


0.4623593 


20 


9 





-30.010574 


-30.174576 


0.8235512 



Table 2. Results for N b = 200 and £ max = 1. 



N = 


= 16 




N = 


= 20 


£ = 


£ = l 


£=0 


£ = \ 


1.9999318 


5.9997504 




1.9999832 


5.9999490 


1.9980654 


5.7710970 




1.9997458 


5.9983572 


0.2281366 


0.0026448 




1.9875134 


2.0084946 


0.0002694 


0.0002720 




0.0045222 


0.0012960 


0.0000120 


0.0000084 




0.0000690 


0.0000552 


0.0000014 


0.0000012 




0.0000058 


0.0000066 


0.0000002 


3.316D-07 




0.0000010 


0.0000002 


6.896D-08 


1.005D-07 




0.0000002 


3.072D-07 



Table 3. Occupation numbers of the HFB minimizer, for Ni = 200 and ^ max = 1- 



5.2.5. Quality of the approximation in terms of the number Nf, of points 

In Table|4] we display the HF and HFB ground state energies for N = 16, £ max = 1 for the 
the optimal value of r max , in terms of the number of discretization points Nf, of the regular 
grid. The convergence is not very fast, but we see that the difference between the HF and 
the HFB energy, as well as the norm of the pairing matrix are of the same order for small 
Nb as they are for larger iV&'s. From this observation we can conclude that it is probably 
not necessary to take Nf, very large in order to decide whether pairing occurs or not. 

5.3. A simplified model for protons and neutrons 

In this section we report on our numerical results concerning a simple model inspired of 
nuclear physics. The interaction between protons and neutrons is not a fundamental law 
of nature because these are composite particles made of quarks, which interact through 
weak, strong and electrostatic forces. A common procedure used in nuclear physics is to 
use empirical forces |41j which involve a small number of parameters which are fitted to 
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N b 


'max 


HF energy 


HFB energy 


difference 


\A\\ 


30 


9 


-19.112314 


-19.117948 


0.005634 


0.425604 


50 


9 


-19.189066 


-19.196012 


0.006946 


0.445728 


100 


9 


-19.222300 


-19.229872 


0.007572 


0.454173 


150 


10 


-19.229494 


-19.237574 


0.008080 


0.461725 


200 


10 


-19.232176 


-19.240308 


0.008132 


0.462363 


250 


10 


-19.233420 


-19.241576 


0.008156 


0.462659 


300 


10 


-19.234094 


-19.242264 


0.008170 


0.462821 


400 


11 


-19.234826 


-19.243068 


0.008242 


0.463905 


500 


11 


-19.235206 


-19.243456 


0.008250 


0.463985 



Table 4. Value of the HF and HFB energies for N = 16 and ^ max = 1 and the (approximate) optimal r mJ1 .. 

experiment or to the known behavior of the model in some limits. The most common forces 
used in practice are the so-called Skyrme [45] and Gogny [20j [21] [16] forces and they depend 
nonlinearly on the state itself. Here we consider an effective force which is fixed and does not 
depend on the quantum state. We also take it spin-independent and isospin-independent. 
Our goal is to test some simple ideas and not to do a real nuclear physics calculation. 

5.3.1. Model 

The nucleon-nucleon potential has been observed to be repulsive at short distances and only 
attractive at medium distances. It decays very fast at infinity. A simple choice to describe 
this is to take 

W(x) = -A - oi e~ bl|2:|2 + a 2 e- b2lxl \ (5.6) 
\x\ 

with a2, ai > 0, bi < b 2 - The constant k is 1 for the proton-proton interaction and 
for the proton-neutron and the neutron-neutron interaction. The other constants usually 
also depend on the isospin (the quantum variable which determines whether a nucleon is a 
neutron or a proton). For simplicity we work here with particles having a definite isospin. 
This means that we assume to have either only protons or only neutrons. In particular we 
want to ask for which strength of the effective force it becomes possible for the protons to 
overcome their Coulomb repulsion and form a bound state. In reality a nucleus is made of 
a certain number of protons and neutrons and one has to use a different HFB state for each 
species. 

In our applications we have chosen for simplicity b\ = 1, &a = 4, a\ = a = 2 02/3. This 
means that the effective force takes the form 

W(x) = ^+a(^e-^ 2 -e-^y (5.7) 

When k = 1, this force is purely repulsive for a ^ 2.87 and it becomes attractive at inter- 
mediate distances for larger a's. The corresponding force is displayed in Figure [7] for a = 1 
and k = 0. 

One can ask several questions concerning the model considered in this section: 
(1) For which value of a does a system of N identical nucleons bind in Hartree-Fock theory? 
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Fig. 7. The effective force — e I 1 + 3e /2 used in our calculation of Section 15.31 The (repulsive) 

Coulomb potential must be added for protons. 

(2) Is there always pairing when there is binding? 

(3) Can pairing effects allow for binding with a smaller a than in HF theory? 

These questions are mostly of academic nature for the very simplified model considered in 
this section. But investigating the same problems with more realistic forces is very important 
from a physical point of view. From a mathematical point of view, nothing seems to be known 
for simple models of the same form as in this section. It is not even known that binding 
always occurs for a large enough with the previous interaction. We hope that our calculations 
will stimulate some further mathematical studies. 

5.3.2. Some computational details 

We always minimize over states having the spin, time-reversal and rotation symmetries. The 
Bach-Frohlich- Jonsson Theorem 12.21 does not apply to the model of this section, hence we 
are making a further approximation here. 

For such symmetric states we have shown in Section 13.3.21 that the energy can be ex- 
pressed in terms of 

p oo />oo 

{ij\mn)t,t> = r 2 dr s 2 ds Xi(r)Xj( r )Xm(s)Xn{s)w e ^(r,s) 
Jo Jo 

where, for the model considered in this section, 
1 f 1 

w t jj (r,s) = -J W (Vr 2 +s 2 -2r^) P/(i) P v (*) dt 
= \ ( Pe(t)Pe>(t)( ; ^ = 

- Ol e -bi(r 2 +s 2 -2rst) + ^ e -6 2 (r 2 +S 2 -2r S t)^ ^ 

For < 1,1' < ^max with 4nax not too large, the Gaussian integrals can be computed exactly 
and it is possible to find the exact expression of we,e>(r, s). 
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The computation of the integral (ij\mn)iji against hat functions is much more tedious, 
however. It is easy to find an exact expression for the Coulomb part, but not so simple for 
the Gaussian part. So we have performed a numerical calculation of these integrals. Since 
we have of the order of (Ni,) 4 integrals, we could not take Nf, too large. The results of the 
previous section indicated that the existence of pairing effects does not depend very much 
on the size of the basis. 

5.3.3. Slow convergence and oscillations of Roothaan 

We have observed that the Roothaan algorithm almost always oscillates, even in the HFB 
case (see some examples in Figures [8j and [TO]) . This is in stark contrast with the results 
of the previous section where the Roothaan algorithm was almost always converging. Some- 
times it very slowly converges in the HF case (see, e.g., Figure |9|). However we have always 
obtained convergence for the HF Roothaan algorithm when a is small enough, that is, when 
it is expected that there is actually no binding. For the case displayed in Figure |H] we have 
a = 20 but the critical a is about ~ 24 (see the next section). 

We conclude that using the ODA is very important for such attractive potentials. The 
same might be true with the more involved forces used in nuclear physics. 

5.3.4. The critical strength 

In finite dimension there is always a minimizer. Saying differently, since the particles are 
trapped in a ball, they always bind. Furthermore we work with rotation-invariant states. 
So, for the true model in infinite dimension, the particles escaping to infinity cannot form a 
bound state of the same kind because they are too far from the (fixed) center of symmetry. 
In this special case they will spread out and have a vanishing energy. 

In Hartree-Fock theory, we conclude that we can detect the loss of binding by looking 
at the last filled HF eigenvalue. When it crosses 0, this corresponds to the last particle 
becoming a scattering state. We can therefore choose as definition for the critical strength a, 
the value at which this eigenvalue is 0. In finite dimensional Hartree-Fock-Bogoliubov theory 
things are less clear and we will not discuss the problem of binding. In our simulations we 
have observed that the HFB ground state density was always rather close to the HF ground 
state density, which suggests that there is binding in HFB as well. 

We have made some calculations for Nb — 50 and N — 4. We found that the critical 
strength is about a c ~ 23.5 in the proton-proton case and a c ~ 17.5 in the neutron-neutron 
case. In Figure [TT] we display the HF and HFB energies as functions of the parameter a, 
for N = 4 and k = 1 (proton- proton case). Figure [T2l is the equivalent result for n = 
(neutron- neutron case). For these calculations we have chosen r max = 3 which is the optimal 
choice for a in a neighborhood of the critical value. Like in the previous section the results 
depend on the radius of the ball in which the system is confined. We see that there is always 
pairing, in the sense that the HFB curve is below the HF curve. This is even more manifest 
in the neutron-neutron case for which the potential is much more attractive than for protons, 
which repel with the Coulomb potential. Also, the norm of the pairing matrix A does not 
vary too much with a, it stays between 0.80 and 0.95 for a in the range 15 ^ a ^ 30, for 
both k = and k = 1. 

From these numerical results we can conclude that pairing seems to happen in this model, 
for any strength a for which there is binding in Hartree-Fock theory. It is an interesting 
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Fig. 9. Same calculation with a = 20 and re = 1 (proton-proton case). The Roothaan algorithm slowly 
converges in the HF case and it oscillates in the HFB case but the two values are very close. 




Fig. 10. Same calculation with a = 20 and re = (neutron- neutron case). The Roothaan algorithm oscillates 
in the HFB case, but the two values are very close. 
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problem to actually prove that pairing always occurs, for instance for a large enough. We 
are not aware of any result of this kind. 




16 18 20 22 24 26 28 30 32 If. 18 20 22 24 28 28 30 32 



Fig. 11. Left: Values of the HF (blue) and HFB (red) ground state energies as functions of a, with JV = 4, 
N b = 50, ^max — 1, fmax — 3 and k — 1 (proton-proton case). The vertical line is the value of a for which 
the last filled eigenvalue vanishes. Right: Values of the two filled HF eigenvalues for the same a. 
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Fig. 12. Same calculations for n = (neutron-neutron case). 
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